!>\file aer_cloud.F
!! This file contains the models of Nenes and Seinfeld (2003) \cite Nenes_2003 ,
!! Fountoukis and Nenes (2005) \cite Fountoukis_2005 and Barahona and
!! Nenes (2008, 2009) \cite Barahona_2008 \cite Barahona_2009 .
!>\author Donifan Barahona  donifan.o.barahona@nasa.gov 

!>\ingroup mg2mg3
!>\defgroup aer_cloud_mod Morrison-Gettelman MP aer_cloud Module
!! according to the models of Nenes & Seinfeld (2003) \cite Nenes_2003, 
!! Fountoukis and Nenes (2005) \cite Fountoukis_2005 
!! and Barahona and Nenes (2008, 2009) \cite Barahona_2008 \cite Barahona_2009 .
!! *** Code Developer: Donifan Barahona donifan.o.barahona@nasa.gov
       MODULE aer_cloud

#ifdef GEOS5
       use MAPL_ConstantsMod, r8 => MAPL_R8
#endif
#ifdef NEMS_GSM
       use physcons, only : MAPL_PI => con_pi
       use machine,  only : r8      => kind_phys
#endif

! according to the models of Nenes & Seinfeld (2003), Fountoukis and Nenes (2005) and Barahona and Nenes (2008, 2009).
! *** Code Developer: Donifan Barahona donifan.o.barahona@nasa.gov
!
!=======================================================================
!

       implicit none
       private

       public :: aerosol_activate
       public :: AerConversion
       public :: AerConversion1
       public :: AerProps
       public :: getINsubset
       public :: init_Aer
       public :: aer_cloud_init
       public :: vertical_vel_variance


       integer, parameter :: nsmx_par=20, npgauss=10



       type :: AerProps
!      sequence
       real, dimension(nsmx_par) :: num, dpg, sig, den, kap
     &,                             fdust, fsoot, forg
       integer :: nmods
       end type AerProps

       interface assignment (=)
       module procedure copy_aer
       end interface


!==================================================================


!
!=================================================================

!
       real*8, dimension(npgauss) :: xgs_par, wgs_par
!

!Global aux variables

       type(AerProps) :: AerPr_base_clean, AerPr_base_polluted

       real*8         :: base_mass_so4_polluted, base_mass_so4_clean,
     &                   base_mass_ss, frac_dust(5), frac_bc, frac_org,
     &                   ahet_dust(5), ahet_bc

!==================================================================

!
!=================================================================

       real*8 ::  sh_ice, doin_ice,vmin_ice, acorr_dust, acorr_bc
     &,           denw_par, cpair_par, dhv_par

       integer :: typeofspec_ice
       logical :: purehet_ice, purehom_ice,  is_gocart

!==================================================================
!
!==========================

       integer, parameter :: maxit_par=100

       real, parameter :: amw_par=18d-3,      ama_par=29d-3
     &,                   grav_par=9.81d0,    rgas_par=8.31d0
     &,                   accom_par=1.0d0,    eps_par=1d-6
     &,                   zero_par=1.0e-20,   great_par=1d20
     &,                   pi_par=mapl_pi,     sq2pi_par=sqrt(pi_par)
!    &,                   pi_par=3.1415927d0, sq2pi_par=sqrt(pi_par)
     &,                   sq2_par=1.41421356237d0
!
     &,                   wmw_ice=018d0,      amw_ice=0.029d0
     &,                   rgas_ice=8.314d0,   grav_ice=9.81d0
     &,                   cpa_ice=1005.1d0,   pi_ice=pi_par
     &,                   depcoef_ice=0.1d0,  thaccom_ice=0.7d0
!
     &,                   To_ice=272.15d0,    Tmin_ice=185.d0
     &,                   Pmin_ice=100.0d0,   Thom=236.0d0
     &,                   rv_ice=rgas_ice/wmw_ice


       CONTAINS

!>\ingroup aer_cloud_mod
!! This subroutine calculates
       subroutine aer_cloud_init()

       real*8  :: daux, sigaux
       integer ::ix

       call AerConversion_base

       acorr_dust = 2.7e7
       acorr_bc   = 8.0e7

       do ix = 1, 5
         daux          = AerPr_base_polluted%dpg(ix)
         sigaux        = AerPr_base_polluted%sig(ix)
         frac_dust(ix) = 0.5d0*(1d0
     &                      - erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))

         ahet_dust(ix) = daux*daux*daux*0.52*acorr_dust
     &                 * exp(4.5*sigaux*sigaux)

       end do


       daux     = AerPr_base_polluted%dpg(12)
       sigaux   = AerPr_base_polluted%sig(12)
       frac_bc  = 0.5d0*(1d0-erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))
       ahet_bc  = daux*daux*daux*0.52*acorr_bc* exp(4.5*sigaux*sigaux)

       daux     = AerPr_base_polluted%dpg(13)
       sigaux   = AerPr_base_polluted%sig(13)
       frac_org = 0.5d0*(1d0-erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))

       end subroutine aer_cloud_init


!>\ingroup aer_cloud_mod
!! This subroutine sets the variables needed for
!!the activation  subroutines and  return the activated droplet and ice number concentration
!>\param tparc_in         T (\f$K\f$)
!!\param pparc_in         P (\f$Pa\f$)
!!\param sigwparc_in      variance of the distribution of updraft velocity (\f$m s^{-1}\f$)
!!\param wparc_ls         mean of the distribution of updraft velocity (\f$m s^{-1}\f$)
!!\param Aer_Props        AerProps structure containing the aerosol properties Aerosol number concentration (\f$Kg^{-1}\f$)
!!\param npre_in          number concentration of prexisting ice crystals (\f$Kg^{-1}\f$)
!!\param dpre_in          mass-weighted diameter of prexisting ice crystals (m)
!!\param ccn_diagr8       array of supersaturations for CCN diagnostics (in-out)
!!\param Ndropr8          Current droplet number concentration (\f$Kg^{-1}\f$)
!!\param cdncr8           Activated cloud droplet number concentration (Kg-1)
!!\param smaxliqr8        Maximum supersaturation w.r.t liquid during droplet activation
!!\param incr8            Nucleated ice crystal concentration (Kg-1)
!!\param smaxicer8        Maximum supersaturation w.r.t. ice during ice nucleation
!!\param nheticer8        Nucleated ice crystal concentration by het freezing (Kg-1)
!!\param INimmr8          Nucleated nc by droplet immersion freezing in mixed-phase clouds (Kg-1)
!!\param dINimmr8         Ice crystal number tendency by immersion freezing (Kg-1 s-1)
!!\param Ncdepr8          Nucleated nc by deposition ice nucleation (Kg-1)
!!\param Ncdhfr8          Nucleated nc by immersion in aerosol (Kg -1)
!!\param sc_icer8         Critical saturation ratio in cirrus
!!\param fdust_immr8      Fraction of immersion mixed-phase ice nuclei that are dust
!!\param fdust_depr8      Fraction of deposition ice nuclei that are dust
!!\param fdust_dhfr8      Fraction of immersion ice nuclei that are dust (not mixed-phase)
!!\param nlimr8           Limiting ice nuclei concentration (m-3)
!!\param use_average_v    .false. integrate over the updraft distribution. True: use the mean vertical velocity
!!\param CCN_param        CCN activation parameterization. 1- Fountoukis and Nenes (2005), 2-Abdul_Razzak and Ghan (2002) (def = 2)
!!\param IN_param         IN activation spectrum (default is 5)
!!\param fd_dust
!!\param fd_soot
!!\param pfrz_inc_r8
!!\param sigma_nuc
!!\param rhi_cell
!!\param nccn
!!
       subroutine aerosol_activate(tparc_in, pparc_in, sigwparc_in,     &
     & wparc_ls, Aer_Props, npre_in, dpre_in, ccn_diagr8, Ndropr8,      &
     & cdncr8, smaxliqr8, incr8, smaxicer8, nheticer8, INimmr8,         &
     & dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8, fdust_immr8, fdust_depr8,  &
     & fdust_dhfr8, nlimr8, use_average_v, CCN_param, IN_param, fd_dust,&
     & fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell,nccn)  
!    & fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell,nccn, lprnt)




       type(AerProps), intent(in) :: Aer_Props

       logical :: use_average_v
!      logical :: use_average_v, lprnt

       real(r8), intent(in) :: tparc_in, pparc_in, sigwparc_in,         &
     & wparc_ls, npre_in, dpre_in, Ndropr8,  fd_soot, fd_dust,          &
     & sigma_nuc, rhi_cell
       integer, intent(in) :: CCN_param, IN_param, nccn

       real(r8), dimension(:), intent(inout) :: ccn_diagr8

       real(r8), intent(out) :: cdncr8, smaxliqr8, incr8, smaxicer8,    &
     & nheticer8, INimmr8, dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8,        &
     & fdust_immr8, fdust_depr8, fdust_dhfr8, nlimr8, pfrz_inc_r8

       type(AerProps) :: Aeraux


       integer :: k, n, I, J, naux


       real*8 :: nact, wparc, tparc,pparc, accom,sigw, smax, antot,     &
     &           ccn_at_s, sigwparc
!      real*8, allocatable, dimension(:) :: smax_diag
       real*8, dimension(nccn) :: smax_diag


       real*8 :: nhet, nice, smaxice, nlim, air_den, frac, norg, nbc,   &
     &           nhom, dorg, dbc, kappa, INimm, dINimm, aux

!Anning Cheng move allocable array here for thread safety, 5/4/2016
!      real*8, allocatable, dimension(:) :: sg_par, tp_par, dpg_par,
       real*8, dimension(nsmx_par)  :: sg_par,    tp_par,    dpg_par
     &,                                sig_par,   vhf_par,   ams_par
     &,                                dens_par,  deni_par,  amfs_par
     &,                                kappa_par, ndust_ice, sigdust_ice
     &,                                ddust_ice
!      real*8, allocatable, dimension(:) :: ndust_ice, sigdust_ice,
!    & ddust_ice
       real*8 :: temp_par,  pres_par
       real*8 :: akoh_par,  alfa_par,  bet2_par
       real*8 aka_par, dv_par, psat_par, dair_par,surt_par,ddry_ice,
     & np_ice,nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
     & g1_ice, g2_ice,gdoin_ice, z_ice, norg_ice, sigorg_ice,
     & dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,Nhet_dep,fdust_imm,fdust_dep,fdust_dhf,
     & waux_ice,fdrop_dust,fdrop_bc,D_preex, N_preex,
     & one_over_tao_preex,P_ice, T_ice,miuv_ice,Nhet_dhf,denice_ice,
     & vpresw_ice,vpresi_ice,denair_ice
       real*8  :: ntot
       integer :: nmodes,act_param,nbindust_ice
       logical :: use_av_v

!=============inputs================
       tparc     = tparc_in
       pparc     = pparc_in
       sigwparc  = sigwparc_in
       miuv_ice  = wparc_ls
       air_den   = pparc*28.8d-3/rgas_par/tparc
       N_preex   = max(npre_in*air_den, zero_par)
       D_preex   = max(dpre_in, 1.0e-9)
       use_av_v  = use_average_v
       act_param = 2
       typeofspec_ice = 5



       smaxicer8   = zero_par
       smaxice     = zero_par
       cdncr8      = zero_par
       smaxliqr8   = zero_par
       incr8       = zero_par
       smaxice     = max(2.349d0-(tparc/259d0) -1.0 , 0.0)
       nheticer8   = zero_par
       nlimr8      = zero_par
       sc_ice      = max(2.349d0-(tparc/259d0), 1.0)
       If (tparc > Thom) sc_ice =1.0

       INimmr8     = zero_par
       dINimmr8    = zero_par
       Ncdepr8     = zero_par
       Ncdhfr8     = zero_par
       fdust_immr8 = zero_par
       fdust_dhfr8 = zero_par
       fdust_depr8 = zero_par
       fdust_imm   = zero_par
       fdust_dhf   = zero_par
       fdust_dep   = zero_par
       pfrz_inc_r8 = zero_par

       nact     = zero_par
       smax     = zero_par
       sc_icer8 = sc_ice

       is_gocart = .true.

       if (sum(Aer_Props%num) <= 1.0e2) then
         return
       end if

       nmodes = max(Aer_Props%nmods, 1)

!      allocate (dpg_par(nmodes))
!      allocate (vhf_par(nmodes))
!      allocate (ams_par(nmodes))
!      allocate (dens_par(nmodes))
!      allocate (sig_par(nmodes))
!      allocate (tp_par(nmodes))
!      allocate (amfs_par(nmodes))
!      allocate (deni_par(nmodes))
!      allocate (sg_par(nmodes))
!      allocate (smax_diag(size(ccn_diagr8)))

!      allocate (kappa_par(nmodes))


       smax_diag = 0.01
       sigw      = zero_par
       do n=1,nmodes
         dpg_par(n)   = zero_par
         vhf_par(n)   = zero_par
         ams_par(n)   = zero_par
         dens_par(n)  = zero_par
         sig_par(n)   = 1d0
         tp_par(n)    = zero_par
         amfs_par(n)  = zero_par
         deni_par(n)  = zero_par
         kappa_par(n) = zero_par
       enddo

       call init_Aer(Aeraux)

!     if (lprnt) write(0,*)' in aero Aer_Props%num='
!    &,Aer_Props%num,' nmodes=',nmodes,' air_den=',air_den
!     if (lprnt) write(0,*)' in aero Aer_Props%kap='
!    &,Aer_Props%kap

       antot = 0.0

       do n=1,nmodes
!        tp_par(n)    = DBLE(Aer_Props%num(n))*air_den
!        dpg_par(n)   = max(DBLE(Aer_Props%dpg(n)), 1.0e-10)
!        sig_par(n)   = DBLE(Aer_Props%sig(n))
!        kappa_par(n) = max(DBLE(Aer_Props%kap(n)), 0.001)
!        dens_par(n)  = DBLE(Aer_Props%den(n))

         tp_par(n)    = Aer_Props%num(n) * air_den
         dpg_par(n)   = max(Aer_Props%dpg(n), 1.0e-10)
         sig_par(n)   = Aer_Props%sig(n)
         kappa_par(n) = max(Aer_Props%kap(n), 0.001)
         dens_par(n)  = Aer_Props%den(n)
         vhf_par(n)   = 3.0
         if (kappa_par(n) > 0.01) then
           ams_par(n) = 18.0e-3*1.7*3.0/kappa_par(n)
         else
           ams_par(n) = 900.0e-3
           tp_par(n)  = 0.0
         endif
         amfs_par(n)  = 1.0
         deni_par(n)  = dens_par(n)
         antot        = antot + tp_par(n)

!     if (lprnt) write(0,*)' n=',n,' tp_par=',tp_par(n),' antot=',antot
!    &,' Aer_Props%num=',Aer_Props%num(n),' kappa_par=',kappa_par(n)
!    &,' air_den=',air_den
       enddo

!      kappa_par = max(kappa_par, 0.001)
!      dpg_par   = max(dpg_par, 1.0e-10)
       temp_par  = max(tparc, 245.0)
       pres_par  = max(pparc, 34000.0)

!      antot = sum(tp_par)
       ntot  = antot

       wparc = max(max(0.8d0*sigwparc, 0.01)+ wparc_ls, 0.01)


       act_param = CCN_param

!============== Calculate cloud droplet number concentration===================

!     if (lprnt) write(0,*)' in aero tparc=',tparc,' antot=',antot
!     if (lprnt) write(0,*)' in aero tp_par=',tp_par(1:nmodes)

       if (tparc > 245.0) then
         if (antot > 1.0) then

           call ccnspec (tparc,pparc,nmodes
     &,      amfs_par,dens_par,deni_par,vhf_par,ams_par
     &,      sg_par,tp_par,akoh_par,surt_par,temp_par,pres_par
     &,      dv_par,act_param,aka_par, psat_par,dair_par
     &,      ntot,dpg_par)

           if (wparc >= 0.005) then
             if (act_param > 1) then

               call arg_activ (wparc,0.d0,nact,smax,nmodes,tp_par
     &,             dpg_par,kappa_par,sig_par,temp_par, pres_par)

             else

               call pdfactiv (wparc,0.d0,nact,smax,nmodes
     &,             alfa_par,bet2_par,akoh_par,sg_par,tp_par
     &,             temp_par, pres_par,aka_par, dv_par, psat_par
     &,             dair_par,ntot,sig_par)

             endif
           endif

           cdncr8    = max(nact/air_den, zero_par)
           smaxliqr8 = max(smax, zero_par)

!     if (lprnt) write(0,*)' in aero cdncr8=',cdncr8,' nact=',nact,
!    &' air_den=',air_den,' wparc=',wparc,' act_param=',act_param
!============ Calculate diagnostic CCN number concentration==================

!          smax_diag = ccn_diagr8

!          do k =1, size (smax_diag)
!          do k =1, nccn
!            call ccn_at_super (smax_diag(k), ccn_at_s,nmodes,
!    &                          sig_par,sg_par,tp_par)
!            ccn_diagr8 (k) = ccn_at_s
!          end do

         end if
       end if


! ==========================================================================================
! ==========================================================================================
!==========================  Ice crystal nucleation parameterization  ======================
! ==========================================================================================
       dbc_ice      = 1.0e-9
       nbc_ice      = zero_par
       norg_ice     = zero_par
       dorg_ice     = 1.0e-9
       sigbc_ice    = zero_par
       sigorg_ice   = zero_par
       ddry_ice     = 1.0e-9
       np_ice       = zero_par
       nin_ice      = 0.
       kdust_ice    = 0.
       kbc_ice      = 0.
       shdust_ice   = 0.
       shbc_ice     = 0.
       effdust_ice  = 0.
       effbc_ice    = 0.
       del1dust_ice = 0.
       si0dust_ice  = 0.
       del1bc_ice   = 0.
       si0bc_ice    = 0.

       naux = 0
       do k = 1, nmodes
         if (kappa_par(k) > 0.1) then
           np_ice   = np_ice   + tp_par(k)
           ddry_ice = ddry_ice + dpg_par(k)
           naux     = naux + 1
         end if
       end do
       ddry_ice = ddry_ice/max(naux , 1)
       frac     = 1.0
       np_ice   = frac*np_ice

!get dust from input structure

       call getINsubset(1, Aer_Props, Aeraux)
       nbindust_ice = max(Aeraux%nmods, 1)

!      allocate(ndust_ice(nbindust_ice))
!      allocate(sigdust_ice(nbindust_ice))
!      allocate(ddust_ice(nbindust_ice))

       ddust_ice = zero_par
       ndust_ice = zero_par
       sigdust_ice = zero_par
       do n=1,nbindust_ice
         ddust_ice(n)   = DBLE(Aeraux%dpg(n))
         ndust_ice(n)   = DBLE(Aeraux%num(n))*air_den
         sigdust_ice(n) = DBLE(Aeraux%sig(n))
       enddo


!Black carbon. Only a single mode considered. Use average size and sigma

       call getINsubset(2, Aer_Props, Aeraux)
       naux      = max(Aeraux%nmods, 1)
       dbc_ice   = DBLE(sum(Aeraux%dpg(1:naux)))/naux
       nbc_ice   = DBLE(sum(Aeraux%num(1:naux)))*air_den
       sigbc_ice = DBLE(sum(Aeraux%sig(1:naux)))/naux



       call getINsubset(3, Aer_Props, Aeraux)
       naux       = max(Aeraux%nmods, 1)
       dorg_ice   = DBLE(sum(Aeraux%dpg(1:naux)))/naux
       norg_ice   = DBLE(sum(Aeraux%num(1:naux)))*air_den
       sigorg_ice = DBLE(sum(Aeraux%sig(1:naux)))/naux


       nhet     = zero_par
       nice     = zero_par
       nlim     = zero_par
       INimm    = zero_par
       dINimm   = zero_par
       Nhet_dep = zero_par
       Nhet_dhf = zero_par
       antot=sum(ndust_ice)+ norg_ice+ nbc_ice+ np_ice

       sigwparc=max(0.01, sigwparc)
       sigwparc=min(5.0, sigwparc)
       waux_ice=max(wparc_ls + sigwparc*0.8, 0.01)


!===========Calculate nucleated crystal number. Follows Barahona and Nenes (2008, 2009)==============

       purehet_ice= .FALSE.
       purehom_ice= .FALSE.


      if (antot > 1.0e2) then
        if (tparc < To_ice) then

          CALL prop_ice(tparc, pparc, denice_ice,ddry_ice
     &,    nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
     &,    g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice
     &,    kdust_ice, kbc_ice, shdust_ice, shbc_ice
     &,    effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice
     &,    si0bc_ice,nbc_ice,D_preex, N_preex
     &,    one_over_tao_preex,P_ice, T_ice,act_param
     &,    ndust_ice,vpresw_ice,vpresi_ice,denair_ice)


          if (tparc > Thom) then

            fdrop_dust = fd_dust
            fdrop_bc   = fd_soot

            if (sum(ndust_ice)+ norg_ice+ nbc_ice .gt. 1.e3) then

              call INimmersion(INimm, dINimm, waux_ice,dbc_ice,sigbc_ice
     &,            nbc_ice,fdust_imm,fdrop_dust,fdrop_bc,ndust_ice
     &,            sigdust_ice,ddust_ice,nbindust_ice
     &,            vpresw_ice,vpresi_ice,T_ice)


              ndust_ice = max(ndust_ice*(1.0-fdrop_dust), 0.0)
              nbc_ice   = max(nbc_ice*(1.0-fdrop_bc), 0.0)


              call IceParam (sigwparc, denice_ice,ddry_ice,np_ice
     &,         nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
     &,         g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
     &,         norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
     &,         kdust_ice, kbc_ice, shdust_ice, shbc_ice
     &,         effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
     &,         del1bc_ice,  si0bc_ice,nbc_ice,one_over_tao_preex
     &,         nhet, nice, smaxice, nlim,Nhet_dep,Nhet_dhf,fdust_dep
     &,         P_ice,T_ice,ndust_ice,sigdust_ice,ddust_ice,nbindust_ice
     &,         use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
            end if

            sc_ice = 1.0

          else

            call IceParam (sigwparc, denice_ice,ddry_ice,np_ice
     &,      nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
     &,      g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
     &,      norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
     &,      kdust_ice, kbc_ice, shdust_ice, shbc_ice
     &,      effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
     &,      del1bc_ice,  si0bc_ice,nbc_ice,one_over_tao_preex
     &,      nhet, nice, smaxice, nlim,Nhet_dep,Nhet_dhf,fdust_dep
     &,      P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
     &,      use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)

          end if

          aux = (sc_ice - rhi_cell)/(sq2_par*sigma_nuc)

          pfrz_inc_r8 = 1.0d0- 0.5d0*(1.0d0+erf(aux))
          pfrz_inc_r8 = min(max(pfrz_inc_r8, 0.0), 0.999)

        end if
      end if


!======================== use sc_ice only for cirrus
       If (tparc > Thom) sc_ice =1.0
!==========================

!All # m-3 except those passed to MG later
       smaxicer8   = min(max(smaxice, zero_par), 2.0)
       nheticer8   = min(max(nhet, zero_par), 1e10)
       incr8       = min(max(nice/air_den, zero_par), 1e10)
       nlimr8      = min(max(nlim, zero_par), 1e10)
       sc_icer8    = min(max(sc_ice, 1.0), 2.0)
       INimmr8     = min(max(INimm, zero_par), 1e10)
       dINimmr8    = min(max(dINimm/air_den, zero_par), 1e10)
       Ncdepr8     = min(max(Nhet_dep, zero_par), 1e10)
       Ncdhfr8     = min(max(Nhet_dhf, zero_par), 1e10)
       fdust_immr8 = min(max(fdust_imm, zero_par), 1e10)
       fdust_depr8 = min(max(fdust_dep, zero_par), 1e10)
       fdust_dhfr8 = min(max(fdust_dhf, zero_par), 1e10)

!      deallocate (ndust_ice)
!      deallocate (sigdust_ice)
!      deallocate (ddust_ice)
!      deallocate (dpg_par)
!      deallocate (vhf_par)
!      deallocate (ams_par)
!      deallocate (dens_par)
!      deallocate (sig_par)
!      deallocate (tp_par)
!      deallocate (amfs_par)
!      deallocate (deni_par)
!      deallocate (sg_par)
!      deallocate (smax_diag)

!      deallocate (kappa_par)



2033    return

       END subroutine aerosol_activate
!


!=======================================================================
!
!>\ingroup aer_cloud_mod
!! This subrotine sets basic properties of the aerosol size distributions when using GOCART aerosol
!!Mass-number conversion based on Barahona at al.(2014) \cite Barahona_2014.
       SUBROUTINE AerConversion_base ()

       integer, parameter    :: NMDM = 20
       real, dimension(NMDM) :: TPI, DPGI, SIGI, DENSI, KAPPAS, FDUST,
     &                          FSOOT, FORG, TPI_aux, DPGI_aux, SIGI_aux

       real:: aux
       integer:: nmod, K

       do k=1,NMDM
         TPI(k)    = 0.0
         DPGI(k)   = 1.0e-9
         SIGI(k)   = 2.0
         DENSI(k)  = 2200.0
         KAPPAS(k) = 0.01
         FDUST(k)  = 0.0
         FSOOT(k)  = 0.0
         FORG(k)   = 0.0
       enddo
       nmod = 13

! Gocart aerosol size distributions for dust

!!!!!!!!!!!!!!======================================
!!!!!!!!!   Dust
!!!!!!!!!!!!!!======================================


       do k=1,5
         SIGI(k)   = log(1.8)
         DENSI(k)  = 1700.0
         KAPPAS(k) = 0.0001
         FDUST(k)  = 1.0
       enddo



       DPGI (1) = 1.46e-6

       DPGI (2) = 2.80e-6

       DPGI (3) = 4.80e-6
!! Dust 4: 3-6
       DPGI (4) = 9.0e-6
!!  Dust 5: 6-10
       DPGI (5) = 16.0e-6

       DO K =1 , 5
         TPI(K) = 6.0/(DENSI(K)*pi_par*exp(4.5*SIGI(K)*SIGI(K))*DPGI(K)*
     &                 DPGI(K)*DPGI(K))
       END DO

!!!!!!!!!!!!!!======================================
!!!!!!!!!   Sea Salt (Using 3 modes based on Barahona et al. GMD. 2014.
!!!!!!!!!!!!!!======================================



       DENSI(6:8)  = 2200.0
       KAPPAS(6:8) = 1.28



       TPI  (6) = 100e6
       DPGI (6) = .01e-6
       SIGI (6) = log(1.6)

       TPI  (7) = 60.0e6
       DPGI (7) = 0.071e-6
       SIGI (7) = log(2.0)

       TPI  (8) = 3.1e6
       DPGI (8) = 0.62e-6
       SIGI (8) = log(2.7)

       aux = 0.
       DO K =6 , 8
         aux = (TPI(K)*DENSI(K)*pi_par*exp(4.5*SIGI(K)*SIGI(K))*DPGI(K)*
     &       DPGI(K)*DPGI(K))/6.0 + aux
       END DO
       base_mass_ss = aux




       KAPPAS(9:11) = 0.65
       DENSI(9:11)  = 1650.0

! Different size distributions for polluted and clean environments


       TPI  (9) = 1.06e11
       DPGI (9) = .014e-6
       SIGI (9) = log(1.8d0)

       TPI  (10) = 3.2e10
       DPGI (10) = 0.054e-6
       SIGI (10) = log(2.16)

       TPI  (11) = 5.4e6
       DPGI (11) = 0.86e-6
       SIGI (11) = log(2.21)

       aux = 0.
       DO K =9, 11
         aux = (TPI(K)*DENSI(K)*pi_par*exp(4.5*SIGI(K)*SIGI(K))*DPGI(K)*
     &          DPGI(K)*DPGI(K))/6.0 + aux
       END DO
       base_mass_so4_polluted = aux




       TPI_aux  (9) = 1.0e9
       DPGI_aux (9) = .016e-6
       SIGI_aux (9) = log(1.6d0)

       TPI_aux  (10) = 8.0e8
       DPGI_aux (10) = 0.067e-6
       SIGI_aux (10) = log(2.1)

       TPI_aux  (11) = 2.0e6
       DPGI_aux (11) = 0.93e-6
       SIGI_aux (11) = log(2.2)


       aux = 0.
       DO K =9, 11
         aux =(TPI_aux(K)*DENSI(K)*pi_par*exp(4.5*SIGI_aux(K)*
     &    SIGI_aux(K))*DPGI_aux(K)*DPGI_aux(K)*DPGI_aux(K))/6.0 + aux
       END DO
       base_mass_so4_clean = aux




       DPGI (12)  = 0.0118*2.e-6
       SIGI (12)  = log(2.00)
       DENSI(12)  = 1600.0
       KAPPAS(12) = 0.0001
       FSOOT(12)  = 1.0
       TPI(12)    = 6.0/(DENSI(12)*pi_par*exp(4.5*SIGI(12)*SIGI(12))*
     &                   DPGI(12)*DPGI(12)*DPGI(12))



       DPGI (13)  = 0.0212*2.e-6
       SIGI (13)  = log(2.20)
       DENSI(13)  = 900.0
       KAPPAS(13) = 0.0001
       FORG(13)   = 1.0
       TPI(13)    = 6.0/(DENSI(13)*pi_par*exp(4.5*SIGI(13)*SIGI(13))*
     &                   DPGI(13)*DPGI(13)*DPGI(13))

       call init_Aer(AerPr_base_polluted)
       call init_Aer(AerPr_base_clean)


       do k=1,nmod
         AerPr_base_polluted%num(k)   = TPI(k)
         AerPr_base_polluted%dpg(k)   = DPGI(k)
         AerPr_base_polluted%sig(k)   = SIGI(k)
         AerPr_base_polluted%kap(k)   = KAPPAS(k)
         AerPr_base_polluted%den(k)   = DENSI(k)
         AerPr_base_polluted%fdust(k) = FDUST(k)
         AerPr_base_polluted%fsoot(k) = FSOOT(k)
         AerPr_base_polluted%forg(k)  = FORG(k)
       enddo
       AerPr_base_polluted%nmods = nmod

       AerPr_base_clean           = AerPr_base_polluted
       AerPr_base_clean%num(9:11) = TPI_aux(9:11)
       AerPr_base_clean%dpg(9:11) = DPGI_aux(9:11)
       AerPr_base_clean%sig(9:11) = SIGI_aux(9:11)

       RETURN
!
       END SUBROUTINE AerConversion_base


!=======================================================================
!
!>\ingroup aer_cloud_mod
!! This subroutine sets the properties of the aerosol distributions
!! Mass-number conversion based on Barahona at al. GMD, 2014. \cite Barahona_2014
       SUBROUTINE AerConversion (aer_mass, AerPr, kappa, SULFATE, ORG,  &
     &                           BCARBON, DUST, SEASALT)


       type(AerProps), intent (out), dimension(:,:,:) :: AerPr


       real, dimension(:,:,:,:), intent(in) :: aer_mass
       real, intent (out), dimension(:,:,:) :: kappa, DUST, SULFATE,    &
     &                                         BCARBON, ORG, SEASALT
       real:: aux, densSO4, densORG, k_SO4, k_ORG, k_SS, tot_mass, dens,
     &        kappa_aux, tx1

       integer :: i,j,k,l
       integer :: im, jm, lm
       type(AerProps) :: AeroAux
       real, dimension(size(aer_mass,4)) :: aer_mass_tmp

       im = size(aer_mass,1)
       jm = size(aer_mass,2)
       lm = size(aer_mass,3)

       call init_Aer(AeroAux)

       do k = 1, lm
         do j = 1, jm
           do i = 1, im
             aer_mass_tmp(:) = aer_mass(i,j,k,:)


             tot_mass  = aer_mass_tmp(11) + aer_mass_tmp (15)
     &                 + aer_mass_tmp (14)
             densSO4   = 1700.0
             densORG   = 1600.0
             k_SO4     = 0.65
             k_ORG     = 0.2
             kappa_aux = 0.65



             if (tot_mass > 2.0e-20) then
               tx1 = 1.0 / tot_mass
               dens = (aer_mass_tmp(11)*densSO4
     &              + sum(aer_mass_tmp(14:15))*densORG) * tx1
               kappa_aux = (aer_mass_tmp(11)*k_SO4
     &                   + sum(aer_mass_tmp(14:15))*k_ORG) * tx1
             else
               dens = 1750.0
               kappa_aux = 0.65
             end if


             if (tot_mass > 5.0e-8) then
               AeroAux = AerPr_base_polluted
               AeroAux%num(9:11) = AeroAux%num(9:11)*tot_mass
     &                           / base_mass_so4_polluted
             else
               AeroAux = AerPr_base_clean
               AeroAux%num(9:11) = AeroAux%num(9:11)*tot_mass
     &                           / base_mass_so4_clean
             end if

             AeroAux%kap(9:11) = max(kappa_aux, 0.001)
             AeroAux%den(9:11) = dens
             SULFATE(i,j,k)    = SUM(AeroAux%num(9:11))
             kappa_aux         = kappa_aux*tot_mass




             AeroAux%num(1:5) = AeroAux%num(1:5) *aer_mass_tmp(1:5)
             kappa_aux        = kappa_aux
     &                        + AeroAux%kap(1)*sum(aer_mass_tmp(1:5))
             DUST(i,j,k)      = sum(AeroAux%num(1:5))


             tot_mass         = sum(aer_mass_tmp(6:10))
             AeroAux%num(6:8) = AeroAux%num(6:8)*tot_mass/base_mass_ss
             kappa_aux        = kappa_aux + AeroAux%kap(6)*tot_mass
             SEASALT(i,j,k )  = sum(AeroAux%num(6:8))


             AeroAux%num(12) = AeroAux%num(12) *aer_mass_tmp(13)
             kappa_aux       = kappa_aux
     &                       + AeroAux%kap(12)*aer_mass_tmp(13)
             BCARBON(i,j,k)  = AeroAux%num(12)


             AeroAux%num(13) = AeroAux%num(13) *aer_mass_tmp(15)
             ORG(i, j, k)    = AeroAux%num(13)
             tot_mass        = sum(aer_mass_tmp)

             if (tot_mass > 0.0) then
               kappa(i,j,k) = kappa_aux/tot_mass
             end if

             AerPr(i,j,k) = AeroAux

           end do
         end do
       end do

       RETURN
!
       END SUBROUTINE AerConversion

!>\ingroup aer_cloud_mod
!! This subroutine sets the properties of the aerosol distributions
       SUBROUTINE AerConversion1 (aer_mass, AerPr)


       type(AerProps),               dimension(:,:) :: AerPr
!      type(AerProps), intent (out), dimension(:,:) :: AerPr


       real, dimension(:,:,:)             :: aer_mass
!      real, dimension(:,:,:), intent(in) :: aer_mass
       real:: aux, densSO4, densORG, k_SO4, k_ORG, k_SS, tot_mass, dens,
     &        kappa_aux, tx1, tx2

       integer :: i,k,l,im,lm
       type(AerProps) :: AeroAux
       real, dimension(size(aer_mass,3)) :: aer_mass_tmp

       im = size(aer_mass,1)
       lm = size(aer_mass,2)

       call init_Aer(AeroAux)

       do k = 1, lm
         do i = 1, im
           aer_mass_tmp(:) = aer_mass(i,k,:)


           tot_mass  = aer_mass_tmp(11) + aer_mass_tmp (15)
     &               + aer_mass_tmp (14)
           densSO4   = 1700.0
           densORG   = 1600.0
           k_SO4     = 0.65
           k_ORG     = 0.2
           kappa_aux = 0.65



           if (tot_mass > 2.0e-20) then
             tx1       = 1.0 / tot_mass
             tx2       = aer_mass_tmp(14) + aer_mass_tmp(15)
             dens      = (aer_mass_tmp(11)*densSO4 + tx2*densORG) * tx1
             kappa_aux = (aer_mass_tmp(11)*k_SO4   + tx2*k_ORG) * tx1
           else
             dens      = 1750.0
             kappa_aux = 0.65
           end if


           if (tot_mass > 5.0e-8) then
             AeroAux = AerPr_base_polluted
             AeroAux%num(9:11) = AeroAux%num(9:11)*tot_mass
     &                         / base_mass_so4_polluted
           else
             AeroAux = AerPr_base_clean
             AeroAux%num(9:11) = AeroAux%num(9:11)*tot_mass
     &                         / base_mass_so4_clean
           end if

           AeroAux%kap(9:11) = max(kappa_aux, 0.001)
           AeroAux%den(9:11) = dens
           kappa_aux         = kappa_aux*tot_mass




           AeroAux%num(1:5) = AeroAux%num(1:5) *aer_mass_tmp(1:5)
           kappa_aux        = kappa_aux
     &                      + AeroAux%kap(1)*sum(aer_mass_tmp(1:5))


           tot_mass         = sum(aer_mass_tmp(6:10))
           AeroAux%num(6:8) = AeroAux%num(6:8)*tot_mass/base_mass_ss
           kappa_aux        = kappa_aux + AeroAux%kap(6)*tot_mass


           AeroAux%num(12) = AeroAux%num(12)*aer_mass_tmp(13)
           kappa_aux       = kappa_aux+AeroAux%kap(12)*aer_mass_tmp(13)


           AeroAux%num(13) = AeroAux%num(13) *aer_mass_tmp(15)
           tot_mass        = sum(aer_mass_tmp)


           AerPr(i,k) = AeroAux

         end do
       end do

       RETURN
!
       END SUBROUTINE AerConversion1

!=======================================================================
!=======================================================================
!=======================================================================

!>\ingroup aer_cloud_mod
!! This subroutine calculates subgrid scale distribution of vertical velocity.
      subroutine vertical_vel_variance(omeg, lc_turb, tm_gw, pm_gw,
     & rad_cool, uwind_gw, tausurf_gw, nm_gw, LCCIRRUS, Nct, Wct, ksa1,
     & fcn, KH, FRLAND, ZPBL, Z, maxkhpbl, wparc_ls, wparc_gw,
     & wparc_cgw, wparc_turb, LTS)


       real(r8), intent(in) :: omeg, tm_gw, lc_turb, rad_cool, uwind_gw,
     & pm_gw
       real , intent(in) :: LCCIRRUS, KH, ZPBL, Z, FRLAND, nm_gw,
     & tausurf_gw, ksa1, fcn, maxkhpbl, Nct, Wct, LTS

       real(r8), intent(out) :: wparc_ls, wparc_gw, wparc_cgw,
     & wparc_turb

       real(r8) :: rho_gw, k_gw, h_gw, c2_gw, dummyW, maxkh, Wbreak



!!!:========= mean V Large scale and radiative cooling
       rho_gw = pm_gw*28.8d-3/rgas_par/tm_gw


       wparc_ls =-omeg/rho_gw/grav_ice + cpa_ice*rad_cool/grav_ice

!!!======== Orographic Gravity gwave (and brackground) initiated (According to Barahona et al. 2013 GMD)

       wparc_gw = 0.0
       k_gw = 2d0*pi_par/LCCIRRUS

       h_gw= k_gw*rho_gw*uwind_gw*nm_gw

       if (h_gw .gt. 0.0) then
       h_gw=sqrt(2.0*tausurf_gw/h_gw)
       else
       h_gw = 0.0
       end if

       Wbreak = 0.133*k_gw*uwind_gw/nm_gw

       wparc_gw = k_gw*uwind_gw*h_gw*0.133
       wparc_gw = min(wparc_gw, Wbreak)
       wparc_gw = wparc_gw*wparc_gw


!!!======== Subgrid variability from Convective Sources According to Barahona et al. 2014 in prep

       wparc_cgw = 0.0
       c2_gw     = (nm_gw+Nct)/Nct
       wparc_cgw = sqrt(ksa1)*fcn*c2_gw*Wct*0.6334
       wparc_cgw = min(wparc_cgw, Wbreak)
       wparc_cgw = wparc_cgw*wparc_cgw

!!!:=========Subgrid scale variance from turbulence


       wparc_turb=KH/lc_turb


       if (LTS > 20.0) then
         wparc_turb=maxkhpbl/lc_turb
       end if


       wparc_turb= wparc_turb*wparc_turb



      end subroutine vertical_vel_variance
!=======================================================================
!=======================================================================
!=======================================================================

!=======================================================================
!>\ingroup aer_cloud_mod
!! This subroutine extracts aerosol props with INactive = typ
      subroutine getINsubset(typ, aerin, aerout)

! typ  : type of aerosol needed: 1 dust, 2 soot, 3 organics
! nbins: number of modes with num>0

       type(AerProps), intent (in)    :: aerin
       type(AerProps), intent (inout) :: aerout
       integer, intent(in)            :: typ

       integer:: k, bin

       call init_Aer(aerout)
       bin = 0

       do k=1, aerin%nmods

         if (typ == 1) then
           if (aerin%fdust(k) > 0.9) then
             bin = bin + 1
             call copy_mode(aerout,aerin, k,bin)
           end if
         elseif (typ == 2) then
           if (aerin%fsoot(k) > 0.9) then
             bin = bin + 1
             call copy_mode(aerout,aerin, k,bin)
           end if
         elseif (typ == 3) then
           if (aerin%forg(k) > 0.9) then
             bin = bin + 1
             call copy_mode(aerout,aerin, k,bin)
           end if
         end if

       end do

       aerout%nmods = max(bin, 1)

      end subroutine getINsubset

!========================subroutines to handle aer strucuture=====================================

!>\ingroup aer_cloud_mod
!! This subroutine handles aer structure.
       subroutine copy_Aer(a,b)

       type (AerProps), intent(out) :: a
       type (AerProps), intent(in)  :: b

       a%num   = b%num
       a%sig   = b%sig
       a%dpg   = b%dpg
       a%kap   = b%kap
       a%den   = b%den
       a%fdust = b%fdust
       a%fsoot = b%fsoot
       a%forg  = b%forg
       a%nmods = b%nmods

       end subroutine copy_Aer

!>\ingroup aer_cloud_mod
!! This subroutine
       subroutine copy_mode(a_out,a_in, mode_in, mode_out)
       type (AerProps), intent(out) :: a_out
       type (AerProps), intent(in) :: a_in
       integer, intent (in) :: mode_in, mode_out

       a_out%num(mode_out)= a_in%num(mode_in)
       a_out%sig(mode_out) = a_in%sig(mode_in)
       a_out%dpg(mode_out) = a_in%dpg(mode_in)
       a_out%kap(mode_out) = a_in%kap(mode_in)
       a_out%den(mode_out) = a_in%den(mode_in)
       a_out%fdust(mode_out) = a_in%fdust(mode_in)
       a_out%fsoot(mode_out) = a_in%fsoot(mode_in)
       a_out%forg(mode_out) = a_in%forg(mode_in)

       end subroutine copy_mode

!>\ingroup aer_cloud_mod
!! This subroutine initialize aerosol properties in MG sheme.
       subroutine init_Aer(aerout)

       type (AerProps), intent(inout) :: aerout
       integer n

       do n=1,nsmx_par
         aerout%num(n)   =  0.
         aerout%dpg(n)   = 1.0e-9
         aerout%sig(n)   = 2.0
         aerout%kap(n)   = 0.2
         aerout%den(n)   = 2200.0
         aerout%fdust(n) = 0.0
         aerout%fsoot(n) = 0.0
         aerout%forg(n)  = 0.0
       enddo
       aerout%nmods = 1

       end subroutine init_Aer



!>\ingroup aer_cloud_mod
!! This subroutine finds the activated droplet number following Abdul-Razzak and Ghan (2000)
!! \cite Abdul_Razzak_2000 .
!!\author  Donifan Barahona
!!donifan.o.barahona@nasa.gov
       subroutine arg_activ (wparc,sigw,nact,smax,nmodes,tp_par,        &
     &  dpg_par,kappa_par,sig_par,temp_par, pres_par)


       real*8, intent(in) :: wparc, sigw
       integer, intent(in):: nmodes
       real*8, intent(out) :: nact, smax
       real*8 :: SMI(nmodes),tp_par(nmodes)
       real*8, dimension(nmodes)::dpg_par,kappa_par,sig_par

       real*8 :: alfa, beta, Akoh, G, T, fi, gi, nui, citai, ui, aux1,
     & PACT, auxx, aux,temp_par, pres_par
       integer :: I, INDEX

       T = min(max(temp_par, 243.0), 323.0)
       alfa=2.8915E-08*(T*T) - 2.1328E-05*T + 4.2523E-03
       beta=exp(3.49996E-04*T*T - 2.27938E-01*T + 4.20901E+01)
       G=exp(-2.94362E-06*T*T*T + 2.77941E-03*T*T - 8.92889E-01*T +
     & 1.18787E+02)
       Akoh= 0.66e-6/T


       DO I = 1, nmodes
       aux =0.667*Akoh/dpg_par(I)
       SMI (I) = (aux*sqrt(aux))/SQRT(2.0*kappa_par(I))
       END DO



       SMI=MAX(SMI, 1.0e-5)
       aux =alfa*wparc*G
       aux = aux*sqrt(aux)/(2.d0*pi_par*980.d0*beta)
       citai = 0.667*Akoh*SQRT(alfa*wparc*G)

       auxx=0.0
       DO INDEX =1, nmodes
       if (tp_par(INDEX) .gt. 0.0) then
       fi=0.5*exp(2.5*sig_par(INDEX))
       gi=1.0+0.25*sig_par(INDEX)
       nui=aux/tp_par(INDEX)
       aux1=fi*((citai/nui)*sqrt(citai/nui)) + gi*(SMI(INDEX)*
     &SMI(INDEX) /(nui+(3.0*citai)))**0.75
       aux1=aux1/(SMI(INDEX)*SMI(INDEX))
       auxx=auxx+aux1

       end if
       end do


       if (auxx .gt. 0.0) then
         smax = 1/sqrt(auxx)
       else
         nact = 0.0
         return
       end if


       auxx = 0.0


       DO INDEX = 1, nmodes
         ui   = sq2_par*log(SMI(INDEX)/smax)/3.0
         aux1 = 0.5*tp_par(INDEX)*(1.0-ERFAPP(ui))
         auxx = auxx + aux1
       END DO

       nact = auxx


       End Subroutine arg_activ

!===================================================================================================

!===================================================================================================




!=====================================================================================================
!>\ingroup aer_cloud_mod
       subroutine ccn_at_super (super,ccn_at_s,nmodes,                  &
     &                          sig_par,sg_par,tp_par)

       integer :: j, I,nmodes
       real*8  :: dlgsg, dlgsp, orism5, ndl, nds, super, ccn_at_s
       real*8, dimension(nmodes) :: sig_par,sg_par,tp_par

       ndl = 0d0

       do j=1, nmodes

         dlgsg = sig_par(j)

         if (sg_par(j) .gt. 0.0) then
           if (super .gt. 0.0) then
             dlgsp = dlog(sg_par(j)/super)
           else
             dlgsp = dlog(sg_par(j)/0.01)
           end if
         else
           dlgsp = 0.0

         end if
         orism5 = 2.d0*dlgsp/(3d0*sq2_par*dlgsg)
         ndl    = (tp_par(j)/2.0)*(1.0-erf(orism5))+ndl

       end do

       ccn_at_s = ndl

       end subroutine ccn_at_super

!=======================================================================

!    subroutine ccnspec
!------------------------------
!     DESCRIPTION
!
! *** subroutine ccnspec
! *** this subroutine calculates the ccn spectrum of the aerosol using
!     the appropriate form of kohler theory and assuming a lognormal aerosol size dist
!
! *** written by athanasios nenes
!!
!      Code Developer
!      Donifan Barahona
!      donifanb@umbc.edu
!=======================================================================
!
       subroutine ccnspec (tparc,pparc,nmodes,
     & amfs_par,dens_par,deni_par,vhf_par,ams_par,
     & sg_par,tp_par,akoh_par,surt_par,temp_par,pres_par
     & ,dv_par,act_param,aka_par,  psat_par,dair_par,
     & ntot,dpg_par)
!

       integer :: nmodes,i,j,k
       real*8  :: tparc, pparc, amfi,denp,vlfs,par1, par2
       real*8, dimension(nmodes)::amfs_par,dens_par,deni_par,
     &                   vhf_par,ams_par,sg_par,tp_par,dpg_par
       real*8 akoh_par,surt_par,temp_par,pres_par,
     &        aka_par, dv_par, psat_par,dair_par,ntot
       integer act_param


       ntot     = zero_par
       temp_par = max(tparc, 245.0)
       pres_par = max(pparc, 34000.0)
!
       call props(pres_par,temp_par,surt_par,dv_par,act_param,
     & aka_par, psat_par,dair_par)


!
       akoh_par = 4d0*amw_par*surt_par/rgas_par/temp_par/denw_par
!

       do k=1,nmodes


       amfi = max(1.0d0-amfs_par(k),zero_par)
       denp = amfs_par(k)*dens_par(k) + amfi*deni_par(k)
       vlfs = amfs_par(k)/dens_par(k)/(amfs_par(k)/dens_par(k)+ amfi/
     &deni_par(k))
       par1 = 4d0*denw_par*ams_par(k)/27d0/vhf_par(k)/denp/amw_par/
     & dpg_par(k)**3


       par1 = par1/vlfs
       par2 = sqrt(max(par1*akoh_par*akoh_par*akoh_par, zero_par))
       sg_par(k)= max(exp(par2) - 1d0, zero_par)
       ntot = ntot + tp_par(k)
       enddo
!


! *** end of subroutine ccnspec ****************************************
!
       return
       end subroutine ccnspec


!=======================================================================
!    -------------------------------
!     DESCRIPTION
!
!>\ingroup aer_cloud_mod
!! This subroutine calculates the ccn activation fraction according
!! to the nenes and seinfeld (2003) parameterization, with
!!     modification for non-contunuum effects as proposed by fountoukis
!!     and nenes (2005). this routine calculates for a pdf of
!!     updraft velocities.
!!
!!\author Athanasios Nenes
!!
!!      Code Developer
!!      Donifan Barahona
!!      donifanb@umbc.edu
       subroutine pdfactiv (wparc,sigw, nact,smax,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,
     &  temp_par, pres_par,aka_par, dv_par, psat_par,
     &  dair_par,ntot,sig_par)
!
!
       integer :: i, isec,nmodes
       real*8 :: tpart, nact, nacti,wparc,smax,ntot, tx1
       real*8 :: pdf, dpnmx,sigw,plimt,probi,whi,wlo, scal, wpi,smaxi,
     & alfa_par,bet2_par,akoh_par,sg_par(nmodes),tp_par(nmodes),
     & tparc,pparc,temp_par, pres_par,aka_par, dv_par, psat_par,
     & dair_par,sig_par(nmodes)
       real*8 :: wpdbg(npgauss), pddbg(npgauss), nadbg(npgauss),
     &           smdbg(npgauss)
!
! *** case where updraft is very small
!
       if (wparc <= 1d-6) then
         smax = 0d0
         nact = 0d0
         isec = 1
         dpnmx = great_par
         return
       endif
!
! *** single updraft case
!
       if (sigw < 1e-10) then
         call activate (wparc,nact,smax,nmodes,
     &         alfa_par,bet2_par,akoh_par,sg_par,tp_par,
     &         temp_par, pres_par,aka_par, dv_par, psat_par,
     &         dair_par,ntot,sig_par)
         wpdbg(1) = wparc
         pddbg(1) = 1.0
         nadbg(1) = nact
         smdbg(1) = smax
!
! *** pdf of updrafts
!
       else
         nact  = zero_par
         smax  = zero_par
         plimt = 1e-3
         probi = sqrt(-2.0*log(plimt*sigw*sq2pi_par))
         whi   = wparc + sigw*probi
         wlo   = 0.05
         scal  = 0.5*(whi-wlo)
         do i=1,npgauss
           wpi = wlo + scal*(1.0-xgs_par(i))
           call activate (wpi,nacti,smaxi,nmodes,
     &          alfa_par,bet2_par,akoh_par,sg_par,tp_par,
     &          temp_par, pres_par,aka_par, dv_par, psat_par,
     &          dair_par,ntot,sig_par)
           tx1  = (wpi-wparc)/ sigw
           pdf  = (1.0/sq2pi_par/sigw) * exp(-0.5*tx1*tx1)
           nact = nact + wgs_par(i)*(pdf*nacti)
           smax = smax + wgs_par(i)*(pdf*smaxi)
           wpdbg(i) = wpi
           pddbg(i) = pdf
           nadbg(i) = nacti
           smdbg(i) = smaxi
           if (pdf < plimt) exit
         enddo
         nact = nact*scal
         smax = smax*scal
       endif
!
       return
!
! *** end of subroutine pdfactiv ****************************************
!
       end subroutine pdfactiv



!=======================================================================

!     DESCRIPTION
!
! *** subroutine activate
! *** this subroutine calculates the ccn activation fraction according
!     to the nenes and seinfeld (2003) parameterization, with
!     modification for non-contunuum effects as proposed by fountoukis
!     and nenes (2005).
!

!=======================================================================
!
       subroutine activate (wparc,ndroplet,smax,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,temp_par, pres_par,
     &  aka_par, dv_par, psat_par,dair_par,ntot,sig_par)

!
       integer :: i,niter,nmodes

       real*8 :: ndrpl, wparc, beta,cf1, cf2,x1,sinteg1,sinteg2,y1, x2,
     &y2,x3,y3, sign,smax, ent_par, ndroplet, nrdpl,bet1_par,
     &  alfa_par,bet2_par,akoh_par,sg_par(nmodes),tp_par(nmodes),
     &  temp_par, pres_par,aka_par, dv_par, psat_par,dair_par,ntot,
     &  wparcel,sig_par(nmodes)
!
! *** setup common block variables
!
       wparcel = wparc

       sinteg1=zero_par
       sinteg2=zero_par
       nrdpl = zero_par

! *** setup constants
!
       alfa_par = grav_par*amw_par*dhv_par/cpair_par/rgas_par/temp_par /
     &temp_par -grav_par*ama_par/rgas_par/temp_par

       bet1_par = pres_par*ama_par/psat_par/amw_par + amw_par*dhv_par*
     & dhv_par/cpair_par/rgas_par/temp_par/temp_par
       bet2_par = rgas_par*temp_par*denw_par/psat_par/dv_par/amw_par/
     & 4d0 +dhv_par*denw_par/4d0/aka_par/temp_par*(dhv_par* amw_par/
     &rgas_par/temp_par - 1d0)
       beta = 0.5d0*pi_par*bet1_par*denw_par/bet2_par/alfa_par/ wparc/
     &dair_par


       cf1 = 0.5d0*(((1.d0/bet2_par)/(alfa_par*wparc))**0.5d0)
       cf2 = akoh_par/3d0

!
! *** INITIAL VALUES FOR BISECTION **************************************

! *** initial values for bisection **************************************
!
       x1 = 1.0d-5

       if (ntot .gt. zero_par) then
       call sintegral (x1,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
       end if

       y1 = (sinteg1*cf1+sinteg2*cf2)*beta*x1 - 1d0
!
       x2 = 1d0
       if (ntot .gt. zero_par) then
       call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
       end if

       y2 = (sinteg1*cf1+sinteg2*cf2)*beta*x2 - 1d0
!
! *** perform bisection *************************************************
!
20     do 30 i=1,maxit_par
       x3 = 0.5*(x1+x2)
!
       if (ntot .gt. zero_par) then
       call sintegral (x3,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
       end if

       y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0

       if (sign(1.d0,y1)*sign(1.d0,y3) .le. zero_par) then
!                                          ! (y1*y3 .le. zero)
       y2 = y3
       x2 = x3
       else
       y1 = y3
       x1 = x3
       endif
!
       if (abs(x2-x1) .le. eps_par*x1) goto 40
       niter = i


30     continue
!
! *** converged ; return ************************************************
!
40      x3 = 0.5*(x1+x2)
!
       if (ntot .gt. zero_par) then
       call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
       end if


       y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0
       smax = x3

       ndroplet=ndrpl
       return
!
! *** end of subroutine activate ****************************************
!
       end subroutine activate


!=======================================================================
!
! *** subroutine sintegral
! *** this subroutine calculates the condensation integrals, according
!     to the population splitting algorithm of nenes and seinfeld (2003)
!     modal formulation according to fountoukis and nenes (2004)
!
! *** written by athanasios nenes
!
!!     Code Developer
!      Donifan Barahona
!      donifanb@umbc.edu
!=======================================================================
!
       subroutine sintegral (spar, summa, sum, summat,wparcel,nmodes,
     &  alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)

!
       integer ::j,i,k,nmodes
       real*8 :: sum, summat, summa, nd(nsmx_par), integ1(nsmx_par),
     & integ2(nsmx_par),alfa_par,bet2_par,akoh_par,sg_par(nmodes),
     & tp_par(nmodes),sig_par(nmodes)
       real*8 :: descr,spar,ratio, ssplt2, ssplt1, sqrt, ssplt, sqtwo,
     & dlgsg,dlgsp, ekth,wparcel
       real*8 :: orism1, orism2, orism3, orism4, orism5
!
! *** here is where the criterion with the descriminant is put. when it
!     is < 0, then set crit2 = .true. otherwise, set the two values of
!     ssplt and continue.
!
       descr = 1d0 - (16d0/9d0)*alfa_par*wparcel*bet2_par* (akoh_par/
     &spar**2)**2
!
       if (descr.le.0d0) then
       ratio = (2.0d7/3.0)*akoh_par*spar**(-0.3824)
       if (ratio.gt.1.0) ratio = 1.0
       ssplt2 = spar*ratio
       else
       ssplt1 = 0.5d0*(1d0-sqrt(descr))
       ssplt2 = 0.5d0*(1d0+sqrt(descr))
       ssplt1 = sqrt(ssplt1)*spar
       ssplt2 = sqrt(ssplt2)*spar
       endif
!
       ssplt = ssplt2
!
! *** calculate integrals
!
       sum = 0
       summat = 0
       summa = 0
!
       sqtwo = sqrt(2d0)

!
       do 999 j = 1, nmodes
!
       dlgsg = sig_par(j)
       dlgsp = dlog(sg_par(j)/spar)
       orism1 = 2.d0*dlog(sg_par(j)/ssplt2)/(3.d0*sqtwo*dlgsg )
       orism2 = orism1 - 3.d0*dlgsg/(2.d0*sqtwo)
       orism3 = 2.d0*dlgsp/(3.d0*sqtwo*dlgsg)-3.d0*dlgsg/(2.d0*sqtwo)
       orism4 = orism1 + 3.d0*dlgsg/sqtwo
       orism5 = 2.d0*dlgsp/(3*sqtwo*dlgsg)
       ekth = exp(9d0/2d0*dlgsg*dlgsg)
       integ1(j) = tp_par(j)*spar*((1-erf(orism1)) - 0.5d0*((sg_par(j)/
     &spar)**2)*ekth*(1-erf(orism4)))

       integ2(j) = (exp(9d0/8d0*dlgsg*dlgsg)*tp_par(j)/sg_par(j))*
     & (erf(orism2) - erf(orism3))
!
! *** calculate number of drops
!
       nd(j) = (tp_par(j)/2.0)*(1.0-erf(orism5))
!
       sum = sum + integ1(j)
       summat = summat + integ2(j)
       summa = summa + nd(j)
999     continue
!
       return
       end subroutine sintegral

!=======================================================================

!     DESCRIPTION
!
! *** subroutine props
! *** this subroutine calculates the thermophysical properties for the CCN activ param
!
! *** written by athanasios nenes
!      Code Developer
!      Donifan Barahona
!      donifanb@umbc.edu

!=======================================================================
!
       subroutine props(pres_par,temp_par,surt_par,dv_par,act_param,
     &  aka_par, psat_par,dair_par)
!
!
       real*8 :: presa,dbig,dlow,coef,aka_par, dv_par, psat_par
       real*8 :: pres_par,temp_par,surt_par,dair_par
       integer act_param
!
       denw_par = 1d3
       dhv_par = 2.25d6
       cpair_par = 1.0061d3
       presa = pres_par/1.013d5
       dair_par = pres_par*ama_par/rgas_par/temp_par
       aka_par = (4.39+0.071*temp_par)*1d-3
       surt_par = sft(temp_par)

       if (act_param .le. 1) then
       dv_par = (0.211d0/presa)*(temp_par/273d0)**1.94
       dv_par = dv_par*1d-4
       dbig = 5.0d-6
       dlow = 0.207683*((accom_par)**(-0.33048))
       dlow = dlow*1d-6



       coef = ((2*pi_par*amw_par/(rgas_par*temp_par))**0.5)

       dv_par = (dv_par/(dbig-dlow))*((dbig-dlow)-(2*dv_par/accom_par) *
     &coef*(dlog((dbig+(2*dv_par/accom_par)*coef)/(dlow+ (2*dv_par/
     &accom_par)*coef))))

       psat_par = vpres(temp_par)*(1e5/1.0d3)



       end if
!
       return
!
! *** end of subroutine props *******************************************
!
       end subroutine props



!PHYSICAL PROPERTIES for Nenes CDNC Activation



!=======================================================================
!
! *** function vpres
! *** this function calculates saturated water vapour pressure as a
!     function of temperature. valid for temperatures between -50 and
!     50 c.
!
! ======================== arguments / usage ===========================
!
!  input:
!     [t]
!     real variable.
!     ambient temperature expressed in kelvin.
!
!  output:
!     [vpres]
!     real variable.
!     saturated vapor pressure expressed in mbar.
!
!=======================================================================
!
       real*8 function vpres (t)
!
       integer ::i
       real*8 :: a(0:6), t,ttemp, vp
       data a/6.107799610e+0, 4.436518521e-1, 1.428945805e-2,
     & 2.650648471e-4, 3.031240396e-6, 2.034080948e-8, 6.136820929e-11/
!
! calculate polynomial (without exponentiation).
!
       ttemp = t-273.0d0
       vp = a(6)*ttemp
       do i=5,1,-1
       vp = (vp + a(i))*ttemp
       enddo
       vpres = vp + a(0)
!
! end of function vpres
!
       return
       end function vpres



!=======================================================================
!
! *** function sft
! *** this function calculates water surface tension as a
!     function of temperature. valid for temperatures between -40 and
!     40 c.
!
! ======================== arguments / usage ===========================
!
!  input:
!     [t]
!     real variable.
!     ambient temperature expressed in kelvin.
!
!  output:
!     [sft]
!     real variable.
!     surface tension expressed in j m-2.
!
!=======================================================================
!
       real*8 function sft (t)
!
       implicit none

!
       real*8 :: t,tpars
!
       tpars = t-273.15d0
       sft = 0.0761-1.55e-4*tpars
!
       return
       end function sft


! ***********************************************************************
!
       subroutine gauleg (x,w,n)
!
! calculation of points and weights for n point gauss integration
! ***********************************************************************
!
       integer :: n,m,i,j
       real*8 :: x(n), w(n),xm,xl,z,p1,p2,p3,pp,z1
       real*8, parameter :: eps_par=1.e-6
       real*8, parameter :: x1=-1.0, x2=1.0
!
! calculation
!
       m=(n+1)/2d0
       xm=0.5d0*(x2+x1)
       xl=0.5d0*(x2-x1)
       do 12 i=1,m
       z=cos(pi_par*(i-.25d0)/(n+.5d0))
 1     continue
       p1=1.d0
       p2=0.d0
       do 11 j=1,n
       p3=p2
       p2=p1
       p1=((2.d0*j-1.)*z*p2-(j-1.d0)*p3)/j
 11    continue
       pp=n*(z*p1-p2)/(z*z-1.d0)
       z1=z
       z=z1-p1/pp
       if(abs(z-z1).gt.eps_par)go to 1
       x(i)=xm-xl*z
       x(n+1-i)=xm+xl*z
       w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
       w(n+1-i)=w(i)
 12    continue
       return
       end subroutine gauleg

!C=======================================================================
!C
!C *** REAL FUNCTION erf (overwrites previous versions)
!C *** THIS SUBROUTINE CALCULATES THE ERROR FUNCTION USING A
!C *** POLYNOMIAL APPROXIMATION
!C
!C=======================================================================
!C
       REAL*8 FUNCTION erf(x)
       REAL*8 :: x
       REAL*8 :: AA(4), axx, y
       DATA AA /0.278393d0,0.230389d0,0.000972d0,0.078108d0/

       y = dabs(dble(x))
       axx = 1.d0 + y*(AA(1)+y*(AA(2)+y*(AA(3)+y*AA(4))))
       axx = axx*axx
       axx = axx*axx
       axx = 1.d0 - (1.d0/axx)
       if(x.le.0.) then
       erf = -axx
       else
       erf = axx
       endif
       RETURN
       END FUNCTION


!=======================================================================
!
! *** real function erf
! *** this subroutine calculates the error function
!
! *** obtained from numerical recipies
!
!=======================================================================
!

!
!      real*8  :: x

!      if(x.lt.0.)then
!        erf=-gammp(.5d0,x**2)
!      else
!        erf=gammp(.5d0,x**2)
!      endif
!      return

!      end function erf

!
!=======================================================================
!
       real*8 function gammln(xx)
!
!=======================================================================
!
!
       integer :: j
       real*8 :: cof(6),stp,half,one,fpf,x,tmp,ser,xx
!
       data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, -
     &1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
       data half,one,fpf/0.5d0,1.0d0,5.5d0/
       x=xx-one
       tmp=x+fpf
       tmp=(x+half)*log(tmp)-tmp
       ser=one
       do 11 j=1,6
       x=x+one
       ser=ser+cof(j)/x
 11    continue
       gammln=tmp+log(stp*ser)
       return
       end function gammln


!
!=======================================================================
!
!      real*8 function gammp(a,x)
!
!=======================================================================
!
!     real*8 :: a,x,gln,gamser,gammcf

!     if(x.lt.0.d0.or.a.le.0.d0)pause
!     if(x.lt.a+1.d0)then
!       call gser(gamser,a,x,gln)
!       gammp=gamser
!     else
!       call gcf(gammcf,a,x,gln)
!       gammp=1.d0-gammcf
!     endif
!     return
!     end function gammp


!
!=======================================================================
!
!     subroutine gcf(gammcf,a,x,gln)
!
!=======================================================================
!
!
!     integer            :: n
!     integer, parameter :: itmax=100
!     real*8, parameter  :: eps_par=3.e-7
!       real*8             :: gln, gold,a0,a1,x,b0,b1,fac,an, &
!                          float,ana,a,anf,g,gammcf
!     gln=gammln(a)
!     gold=0.
!     a0=1.
!     a1=x
!     b0=0.
!     b1=1.
!     fac=1.
!     do 11 n=1,itmax
!       an=float(n)
!       ana=an-a
!       a0=(a1+a0*ana)*fac
!       b0=(b1+b0*ana)*fac
!       anf=an*fac
!       a1=x*a0+anf*a1
!       b1=x*b0+anf*b1
!       if(a1.ne.0.)then
!         fac=1./a1
!         g=b1*fac
!         if(abs((g-gold)/g).lt.eps_par)go to 1
!         gold=g
!       endif
!1    continue
!     pause 'a too large, itmax too small'
!     gammcf=exp(-x+a*log(x)-gln)*g
!     return
!     end subroutine gcf


!
!=======================================================================
!sft
!     subroutine gser(gamser,a,x,gln)
!
!=======================================================================
!

!     integer  :: n
!     integer, parameter :: itmax=100
!     real*8, parameter  :: eps_par=3.e-7
!       real*8  :: gln,x,gamser,a,ap,sum,del,abs

!     gln=gammln(a)
!     if(x.le.0.)then
!       if(x.lt.0.)pause
!       gamser=0.
!       return
!     endif
!     ap=a
!     sum=1./a
!     del=sum
!     do 11 n=1,itmax
!       ap=ap+1.
!       del=del*x/ap
!       sum=sum+del
!       if(abs(del).lt.abs(sum)*eps_par)go to 1
!1    continue
!     pause 'a too large, itmax too small'
!     gamser=sum*exp(-x+a*log(x)-gln)
!     return
!     end subroutine gser

! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!

!

!*************************************************************************
!
! ICE NUCLEATION PARAMETERIZATION FILES START HERE
!
! ************************************************************************
!
!
!======================================================================
!
!       Code Developer
!       Donifan Barahona, GA TECH
!       donifan@gatech.edu
!    -------------------------------
!     DESCRIPTION
!

!***********************************************************
!** Parameterization  of ICE crystal number concentration
!** for large scale models.
!** Donifan Barahona, Athanasios Nenes
!   JGR, 111, D11211,  2008
!   ACP, 9, 369-381,   2009
!   ACP  9, 5933-5948, 2009
!   Homogeneoeus and heterogeneous nucleation considered
!*** SI units unless otherwise specified.
!
!
! *** WRITTEN BY DONIFAN BARAHONA
!
!=======================================================================



       SUBROUTINE IceParam (sigma_w, denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
     & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,one_over_tao_preex,
     & nhet, nice, smax, nlim,Nhet_dep,Nhet_dhf,fdust_dep
     & ,P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)

       real*8, intent(in) :: sigma_w,denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
     & g1_ice, g2_ice,gdoin_ice, z_ice, lambda_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,P_ice, T_ice,
     & vpresw_ice,vpresi_ice,denair_ice
       real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice


       real*8, intent(out) :: nice, nhet, smax, nlim,Nhet_dep,Nhet_dhf
     & ,fdust_dep

       real*8 :: wpar_ice, sigmav_ice,vmax_ice,miuv_ice
       integer, intent(in)::nbindust_ice
       logical use_av_v




       vmin_ice=0.005d0
       vmax_ice=2.5d0
       sigmav_ice=sigma_w
       vmax_ice= max(min(miuv_ice+(4d0*sigmav_ice), vmax_ice),
     & vmin_ice +0.01)

       if ((sigmav_ice .lt. 0.05) .or. (T_ice .gt. Thom)) then
       use_av_v= .TRUE.
       end if

       if (vmax_ice .gt. vmin_ice + 0.01) then
       if (use_av_v) then
       wpar_ice=min(max(miuv_ice + sigma_w*0.8, vmin_ice), vmax_ice)
       CALL nice_param(wpar_ice, denice_ice,ddry_ice,np_ice,nin_ice,
     & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
     & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
     & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice,denair_ice)
       else
       call nice_Vdist(denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
     & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
     & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
     & ,vpresw_ice,vpresi_ice,denair_ice)
       end IF
       else
       nice = zero_par
       nhet = zero_par
       nlim = 1.0e10
       smax = zero_par
       end if


       return
       END subroutine IceParam


!*************************************************************
!    Subroutine nice_Vdist. Calculates the ice crystal number concentration
!    at the maximum supersaturation using a PDF of updraft using a
!    sixth order Gauss-Legendre quadrature
!     Inputs:  T, and P all SI units)
!    Output NC, smax, nhet, nlim (m-3)
!    Barahona and Nenes, JGR, D11211 (2008) and ACPD, 15665-15698, (2008)
!   Written by Donifan Barahona
!************************************************************


       subroutine nice_Vdist(denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
     & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
     & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice,denair_ice)

       real*8 :: quadx(6), wpar, sum1, quadw(6), dp, sum2, sum3, sum4,
     & sum5, x1, x2
       real*8, intent(out) :: nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf
     &  ,fdust_dep
       real*8  wpar_icex,denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & normv_ice,sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,sc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,one_over_tao_preex,P_ice, T_ice,
     & vpresw_ice,vpresi_ice,denair_ice
       real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice

       INTEGER :: INDEX,nbindust_ice

       DATA quadx/0.23861918d0, -0.23861918d0, 0.66120939d0, -
     &0.66120939d0, 0.93246951d0, -0.93246951d0/

       DATA quadw/0.46791393d0, 0.46791393d0, 0.36076157d0,
     & 0.36076157d0, 0.17132449d0, 0.17132449d0/



       x1=(vmin_ice-miuv_ice)/(sq2_par*sigmav_ice)
       x2=(vmax_ice-miuv_ice)/(sq2_par*sigmav_ice)

       x2=max(x1 +0.01, x2)
       normv_ice=(ERFAPP(x2)-ERFAPP(x1))*0.5d0

       sum1=0d0
       sum2=0d0
       sum3=0d0
       sum4=0d0
       sum5=0d0

       DO INDEX =1, 6
       wpar=max(0.5d0*(((vmax_ice-vmin_ice)*quadx(INDEX)) +(vmax_ice+
     &      vmin_ice)), 0.01)


       CALL nice_param(wpar, denice_ice,ddry_ice,np_ice,nin_ice,
     & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
     & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
     & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice,denair_ice)

       CALL gausspdf(wpar, dp, sigmav_ice,miuv_ice,normv_ice)
       sum1=sum1+(nice*dp*quadw(INDEX))
       sum2=sum2+(smax*dp*quadw(INDEX))
       sum3=sum3+(nhet*dp*quadw(INDEX))
       sum4=sum4+(nlim*dp*quadw(INDEX))
       sum5=sum5+(sc_ice*dp*quadw(INDEX))



       END DO
       nice=sum1*(vmax_ice-vmin_ice)*0.5d0
       smax=sum2*(vmax_ice-vmin_ice)*0.5d0
       nhet=sum3*(vmax_ice-vmin_ice)*0.5d0
       nlim=sum4*(vmax_ice-vmin_ice)*0.5d0
       sc_ice=sum5*(vmax_ice-vmin_ice)*0.5d0
       RETURN


       END subroutine nice_Vdist

!*************************************************************

!*************************************************************
       real*8 function ERFAPP(x)

       real*8, intent(in) :: x
       real*8 :: a
       a = x*x*(1.27324d0+(0.147d0*x*x))/(1d0+(0.147d0*x*x))
       ERFAPP = SQRT(1d0-exp(-a))

       if (x .lt. 0.0) then
         ERFAPP = - ERFAPP
       end if

       end function ERFAPP



!*************************************************************
!    Subroutine Het_freezing. Use only for mixed phase clouds . Inputs: Wpar, T, and P all SI units)
!    Output Nc=Nhet (m-3)
! !    Written by Donifan Barahona
!************************************************************


!      SUBROUTINE Het_freezing(nhet, nice, smax)

!      real*8, intent(out) :: nice, nhet, smax
!      real*8, intent(in) :: np_ice

!      real*8 :: I, SX, NHET_, DSH

!      nhet=zero_par
!      smax=zero_par
!      nice=zero_par

!      SX=vpresw_ice/vpresi_ice
!      call INSPEC_ice(SX-1.0, NHET_, DSH,np_ice)
!      sc_ice=1.0
!      nhet=max(NHET_, zero_par)
!      nice=nhet
!      smax=max(SX-1.0, zero_par)


!      END SUBROUTINE Het_freezing


!*************************************************************
!    Subroutine nice_param. This is the implementation of the Barahona and Nenes(2008, 2009)
!    parameterization
!    Written by Donifan Barahona
!************************************************************


       subroutine nice_param(wpar_icex,denice_ice,ddry_ice,np_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
     & nice, smax, nhet, nlim_,Nhet_dep,Nhet_dhf,fdust_dep,
     & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice,denair_ice)

       real*8, intent(in) :: wpar_icex,denice_ice,np_ice,nin_ice,
     &  alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
     &  sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
     &  norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice, ddry_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,one_over_tao_preex,P_ice, T_ice,
     & vpresw_ice,vpresi_ice,denair_ice
       real*8, intent(out) :: nice, smax, nhet, nlim_,Nhet_dep,Nhet_dhf
     & ,fdust_dep

       real*8 :: AUX1, AUX2, G, DPMAX, MIU, MONVOL, FDS, NLIM, DLIM,
     & DSTAR, DS, NSTAR, NHOM, FC, PHIDO, AUXNC, SIZECORR, DSH, NHET_,
     & F1, F2, SAUX, SUP, SLOW, SHALF, FHALF, DPMIN, GAM, wpar_ice,
     & preex_effect, swsat, sstep,sc_ice,lambda_ice
       real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice

       integer :: INDEX, maxiter_s,nbindust_ice





       preex_effect=1.0-(one_over_tao_preex*(shom_ice)/alfa_ice/
     & ((shom_ice+1.0))/wpar_icex)
       swsat = vpresw_ice/vpresi_ice -1d0


       if (preex_effect .le. 0.0) then
       nhet=0d0
       NHOM=0d0
       smax=shom_ice
       DSH =0.d0
       FDS=1.d0

       sc_ice = shom_ice + 1.d0
       nice = 0.d0
       nlim_=0d0
       return
       else

       wpar_ice = wpar_icex*preex_effect
       end if




       if (np_ice .gt. 1.0) then

       MONVOL=np_ice*1.0d-6*ddry_ice*ddry_ice*ddry_ice
       AUX1=1.6397d-14*T_ice-3.1769d-12
       DPMAX=AUX1*(MONVOL**(-0.373d0))*(wpar_ice**(-0.05))
       IF (DPMAX.GT.1.0d-4) THEN
       DPMAX=1.0d-4
       END IF

       else
       DPMAX=1.0d-4

       endif


       DPMIN=dliq_ice+(0.02/sqrt(alfa_ice*wpar_ice*g1_ice))
       DPMAX=max(DPMAX,DPMIN)

       AUX1=DPMAX-dliq_ice
       AUX2=dlog((g2_ice+(g1_ice*DPMAX))/(g2_ice+(g1_ice*dliq_ice)))
       G=1.3346d0*((g1_ice*AUX1)-(g2_ice*AUX2))/(AUX1*g1_ice*g1_ice)
       lambda_ice=lambda_ice/sqrt(wpar_ice)
       AUX1 = g1_ice*alfa_ice*wpar_ice
       NSTAR=(AUX1*SQRT(AUX1))/beta_ice/z_ice/sq2_par

       GAM=g2_ice/g1_ice




       FDS=1d0
       NHOM=0d0



       if (typeofspec_ice .ge. 0d0) then

       call INSPEC_ice(shom_ice,NHET_,DSH,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     &   kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
     & ,vpresw_ice,vpresi_ice)
       SIZECORR=EXP(-2d0/lambda_ice/shom_ice)
       DSTAR=((4d0*DSH*DSH/3d0)+(2d0*DSH*(shom_ice-DSH))) /(shom_ice-
     &DSH+1d0)

       if (DSTAR .gt. 0.0) then
       NLIM=min(NSTAR*(shom_ice+1d0)/shom_ice/sqrt(DSTAR)/SIZECORR,
     & 1.e10)
       else
       NLIM=1.e10
       end if

       else
       DSH=shom_ice-sh_ice
       DSTAR=((4d0*DSH*DSH/3d0)+(2d0*DSH*(shom_ice-DSH))) /(shom_ice-
     &DSH+1d0)
       DLIM=-GAM+sqrt((GAM*GAM)+(2d0*DSTAR/g1_ice/alfa_ice/wpar_ice))
       NLIM=alfa_ice*wpar_ice*(shom_ice+1d0)/z_ice/beta_ice/shom_ice
       NLIM=NLIM*((g1_ice*DLIM)+g2_ice)/DLIM/DLIM
       NHET_=nin_ice
       end if


       nlim_=min(NLIM, 1d10)
       nlim_=max(NLIM, 1d-6)

       if (NHET_ .gt. 0.0) then
       AUX1 =NHET_/nlim_
       FDS=1d0-(AUX1*SQRT(AUX1))
       else
       FDS = 1d0
       end if



       if (purehom_ice) then
       FDS=1.0d0
       end if

       if ((purehet_ice) .or. (T_ice .GE. Thom)) then
       FDS = 0d0
       end if


       IF (FDS .GE. 1.e-10) THEN

       MIU=FDS*alfa_ice*(shom_ice+1d0)*wpar_ice*koft_ice/shom_ice

       PHIDO=sqrt(pi_ice*G/MIU/2d0)*(G/MIU)
       AUXNC=2d0*denair_ice/beta_ice/koft_ice/denice_ice/pi_ice/np_ice
       FC=AUXNC/PHIDO



       if (np_ice .gt.0d0) then

       IF (FC .le. 0.6d0) then
       NHOM=np_ice*EXP(-FC)*(1.0d0-EXP(-FC))
       else
       NHOM=np_ice/(1d0+EXP((9d0-2d0*FC)/7d0))
       end if

       else
       NHOM=0d0
       end if

       smax=shom_ice
       nhet=NHET_
       if (purehom_ice) NHET_ = 0.d0


       ELSE


       NHOM = 0d0


       smax=0d0
       nhet=0d0
       SAUX=0.01d0

       if (typeofspec_ice .lt. 0d0) SAUX=sh_ice+0.00000000001d0

       F1= FINDSMAX(SAUX,DSH,
     & NHET_,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
     &  gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
     &  beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)
       F2=1.0
       sstep = 0.05d0
       maxiter_s = int(1d0/sstep) + 1

       do INDEX =1, maxiter_s

       if (SAUX .ge. swsat) then
       F2=F1
       SAUX = swsat
       exit
       end if


       SAUX=SAUX+sstep
       F2=FINDSMAX(SAUX,DSH,
     & NHET_,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
     &  gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
     &  beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)
       IF (F2*F1 .lt. 0d0) exit

       end do

       if (F2*F1 .gt.0d0) then
       nhet=0d0
       smax=SAUX
       else

       if (SAUX .lt. swsat) then


       SUP=SAUX
       SLOW=SAUX-(sstep + 0.001d0)

       DO INDEX=1,50
       SHALF=0.5d0*(SUP+SLOW)
       FHALF=FINDSMAX(SHALF,DSH,
     & NHET_,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
     &  gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
     &  beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)

       IF (SIGN(1.d0,F1)*SIGN(1.d0,FHALF) .LE. 0d0) THEN
       F2 = FHALF
       SUP = SHALF
       ELSE
       F1 = FHALF
       SLOW = SHALF
       ENDIF


       IF (ABS(SLOW-SUP) .LE. 5d-3) exit
       END DO
       else
       SHALF = swsat
       end if

       smax=SHALF

       end if

       if ((smax .gt. shom_ice) .and.(T_ice .LE. Thom)) smax=shom_ice

       if (typeofspec_ice .ge. 0d0) then
       call INSPEC_ice(smax, nhet,DSH,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice)
       else
       nhet=nin_ice
       end if

       END IF


       nice=NHOM+nhet
       sc_ice=max(smax+1.0-DSH, 1.0)


      if (.true.) then

       if (FDS .gt. 0.0) then
       sc_ice = (shom_ice+1.0)*FDS + sc_ice*(1.0-FDS)
       end if
       end if

       sc_ice=min(shom_ice+1.0, sc_ice)


       return

       END subroutine nice_param
!*************************************************************
       real*8 function FINDSMAX(SX,DSH,
     & NHET_,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
     &  gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
     &  beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)
       real*8, intent(in) :: SX
       real*8 :: tao
       real*8 ::  DSTAR, SIZECORR, DSH, NSTAR,DLIM
       integer nbindust_ice
       real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
       real*8 shom_ice,NHET_,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  lambda_ice,
     &  gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
     &  beta_ice,nin_ice,vpresw_ice,vpresi_ice

       if (typeofspec_ice .ge. 0d0) then

       call INSPEC_ice(SX,NHET_,DSH,np_ice,norg_ice,sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     &  ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice)
       SIZECORR=EXP(-2d0/lambda_ice/SX)
       DSTAR=((4d0*DSH*DSH/3d0)+(2d0*DSH*(SX-DSH)))/(SX-DSH+1d0)
       DSTAR=DSTAR+(gdoin_ice*alfa_ice*wpar_ice)


       tao=NHET_*SIZECORR*SX*sqrt(DSTAR)/(SX+1d0)/NSTAR


       else
       DSH=SX-sh_ice
       DSTAR=((4d0*DSH*DSH/3d0)+(2d0*DSH*(SX-DSH))) /(SX-DSH+1d0)
       DLIM=-GAM+sqrt((GAM*GAM)+(2d0*DSTAR/g1_ice/alfa_ice/wpar_ice))
       tao=alfa_ice*wpar_ice*(SX+1d0)/z_ice/beta_ice/SX
       tao=tao*((g1_ice*DLIM)+g2_ice)/DLIM/DLIM/nin_ice

       end if


       FINDSMAX=1d0-tao

       end function FINDSMAX


!*************************************************************
!    Function VPRESWATER. Calculates the saturated vapor pressure
!    of water (Pa) according to Murphy & Koop (2005)
!    T in K (173.15-373.15)
!************************************************************

       real*8 function VPRESWATER_ice(T)

       real*8, intent(in) :: T
       real*8 :: A(0:9)

       DATA A/54.842763d0, -6763.22d0, -4.21d0, 0.000367d0, 0.0415d0,
     & 218.8d0, 53.878d0, -1331.22d0, -9.44523d0, 0.014025d0/


       VPRESWATER_ice = A(0)+(A(1)/T)+(A(2)*LOG(T))+(A(3)*T)+
     & (TANH(A(4)*(T-A(5)))*((A(6)+(A(7)/T))+ (A(8)*LOG(T))+ (A(9)*T)))

       VPRESWATER_ice=EXP(VPRESWATER_ice)

       return
       END function VPRESWATER_ice

!*************************************************************
!    Function VPRESICE. Calculates the saturated vapor pressure
!    of ice (pa) according to Murphy & Koop (2005)
!    T in K (>110)
!************************************************************

       real*8 function VPRESICE(T)

       real*8, intent(in) :: T
       real*8 :: A(0:3)

       DATA A/9.550426d0, -5723.265d0, 3.53068d0, -0.00728332d0/


       VPRESICE = A(0)+(A(1)/T)+(A(2)*LOG(T))+(A(3)*T)
       VPRESICE=EXP(VPRESICE)

       return
       END function VPRESICE

!*************************************************************
!    Function DHSUB. Calculates the latent heat of sublimation
!    of ice (J/Kg) according to Murphy & Koop (2005)
!    T in K (>30)
!************************************************************

       real*8 function DHSUB_ice(T)

       real*8, intent(in) :: T
       real*8 :: A(0:4)


       DATA A/46782.5d0, 35.8925d0, -0.07414d0, 541.5d0, 123.75d0/

       DHSUB_ice = A(0) + (A(1) * T) + (A(2)*T*T) + (A(3) * EXP(-((T/
     & A(4))**2)))

       DHSUB_ice=1000d0*DHSUB_ice/18d0
       return
       END function DHSUB_ice

!*************************************************************
!    Function ICEDENSITY. Calculates the DENSITY OF ICE
!    of ice (Kg/m3) according to PK97
!    T in K (>30)
!************************************************************

       real*8 function DENSITYICE(T)

       real*8, intent(in) :: T
       real*8 :: A(0:2), TTEMP

       DATA A/0.9167d0, -1.75d-4, -5.0d-7/

       TTEMP=T-273d0

       DENSITYICE= 1000d0*(A(0)+(A(1)*TTEMP)+(A(2)*TTEMP*TTEMP))
       return
       END function DENSITYICE

!*************************************************************
!    Function WATDENSITY. Calculates the DENSITY OF ICE
!    of liquid water (Kg/m3) according to PK97
!    T in K (>240)
!************************************************************

       real*8 function WATDENSITY_ice(T)

       real*8, intent(in) :: T
       real*8 :: A(0:6), TTEMP, WATDENSITY
       INTEGER :: I


       DATA A/0.99986d0, 6.690d-5, -8.486d-6, 1.518d-7, -6.9984d-9, -
     &3.6449d-10, -7.497d-12 /

       TTEMP=T-273D0

       IF (TTEMP .le. -40d0) THEN
       TTEMP=-40d0
       END IF

       WATDENSITY=A(6)*TTEMP

       IF (T .GE. 240.0) THEN
       DO I=5,1, -1
       WATDENSITY= (WATDENSITY+A(I))*(TTEMP)
       ENDDO
       WATDENSITY=WATDENSITY + A(0)
       ELSE
       WATDENSITY=0.979d0
       END IF

       WATDENSITY=WATDENSITY*1000d0
       WATDENSITY_ice=WATDENSITY
       return
       END function WATDENSITY_ice


!*************************************************************
!    Subroutine PROPERTIES. Set physical an thermodynamic
!    properties at T and P for ice param
!************************************************************


       SUBROUTINE prop_ice(T, P, denice_ice,ddry_ice,
     & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
     & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,
     & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,D_preex, N_preex,one_over_tao_preex,
     & P_ice, T_ice,act_param,ndust_ice,vpresw_ice,vpresi_ice,
     & denair_ice)

       real*8, intent(in)  :: T, P,ddry_ice,nbc_ice,D_preex, N_preex
       real*8, intent(out) :: alfa_ice,beta_ice,shom_ice,
     &  koft_ice, dliq_ice,g1_ice, g2_ice,gdoin_ice, z_ice,
     &  lambda_ice,kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,one_over_tao_preex,P_ice, T_ice

       real*8 :: AUX, AUX1, AUX2, SW, fice, mice, Tc, hdust, hbc, b0,
     &           b1, b2, b3, x, T0bc, T0dust, gam, gamma,NIN_ICE,
     &           Tr, vw, den_m, Eact, Toact, Dact, acc, n1, Siw, rgo,
     &           ngo, mw_molec,dhs_ice,denwat_ice,denice_ice,vpresw_ice,
     &           vpresi_ice,denair_ice,diff_ice,aircond_ice

       real*8,dimension(:) :: ndust_ice
       integer, intent(in) :: act_param


       T_ice    = min(max(T, Tmin_ice), To_ice)
       P_ice    = max(P, Pmin_ice)
       dliq_ice = 1.0d-6


!      rv_ice     = rgas_ice/wmw_ice
       dhs_ice    = DHSUB_ice(T_ice)
       vpresw_ice = VPRESWATER_ice(T_ice)
       vpresi_ice = VPRESICE(T_ice)
       denice_ice = DENSITYICE(T_ice)
       denwat_ice = WATDENSITY_ice(T_ice)
       denair_ice = P_ice*amw_ice/rgas_ice/T_ice


       diff_ice=(0.211d0*101325d0/P_ice)*((T_ice/273d0)**1.94d0)*1.0d-4
       AUX1=1.0e-3*(4.39d0+0.071d0*T_ice)


       AUX2=(2d0*AUX1/(thaccom_ice*1.0d-6*denair_ice*cpa_ice)) *((58.0d-
     &3*pi_ice/(rgas_ice*T_ice))**0.5d0)

       aircond_ice=AUX1/(1.d0+AUX2)




       AUX1=grav_ice*dhs_ice/rv_ice/T_ice/T_ice/cpa_ice
       AUX2=grav_ice*amw_ice/rgas_ice/T_ice
       alfa_ice=AUX1-AUX2
       beta_ice=amw_ice*P_ice/wmw_ice/vpresi_ice
       gamma=1.5d0*dhs_ice*dhs_ice/rv_ice/T_ice/T_ice/cpa_ice

       beta_ice=beta_ice+gamma



       shom_ice=2.349d0-(T_ice/259d0)
       SW=shom_ice*vpresi_ice/vpresw_ice
       shom_ice=shom_ice-1d0
       koft_ice=(0.0240d0*T_ice*T_ice)-(8.035d0*T_ice)+934.0d0





       if (SW .lt. 0.99) then
       AUX1=(1d0/(1d0-SW))-1.1764d0
       else
       AUX1=(1d0/0.01)-1.1764d0
       end if
       dliq_ice=ddry_ice*0.9344d0*(AUX1**0.333)



       AUX1=denice_ice*rv_ice*T_ice/vpresi_ice/diff_ice
       AUX2=dhs_ice*denice_ice/aircond_ice/T_ice
       AUX2=AUX2*((dhs_ice/rv_ice/T_ice)-1.0d0)
       g1_ice=(AUX1+AUX2)/4.0d0


       one_over_tao_preex = beta_ice*denice_ice*pi_ice*0.5* D_preex*
     &N_preex/g1_ice/denair_ice



       g2_ice=denice_ice*rv_ice*T_ice/2.0d0/vpresi_ice/depcoef_ice
       g2_ice=g2_ice*((2.0d0*pi_ice/rv_ice/T_ice)**0.5d0)

       doin_ice=1.0d-6
       gdoin_ice=(g1_ice*0.5d0*doin_ice*doin_ice)+(g2_ice*doin_ice)
       z_ice=denice_ice*pi_ice/2.0d0/denair_ice

       gam=g2_ice/g1_ice
       lambda_ice=1d0/sqrt(alfa_ice*g1_ice*gam*gam)





       if (typeofspec_ice .lt. 0) then
       sh_ice=0.3d0
       nin_ice=(sum(ndust_ice)+nbc_ice)*0.05d0
       elseif (typeofspec_ice .eq. 3) then

       shdust_ice = 0.2d0
       effdust_ice=0.6d0
       shbc_ice = 0.35d0
       effbc_ice=0.05d0
       mice = 0.96d0
       fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
       kdust_ice=koft_ice*fice
       mice = 0.76d0
       fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
       kbc_ice=koft_ice*fice
       elseif (typeofspec_ice .eq. 4) then

       Tc=T_ice-273.15d0
       hdust=0.15d0
       T0dust=-40d0
       b0=-1.0261d0; b1=3.1656d-3; b2=5.3938d-4; b3=8.2584d-6
       x=b0+(b1*Tc)+(b2*Tc*Tc)+(b3*Tc*Tc*Tc)
       si0dust_ice=1d0+(10d0**x)
       del1dust_ice=cubicint_ice(Tc, T0dust, T0dust+5d0, 1d0, hdust)
       hbc=0d0
       T0bc=-50d0
       b0=0.5652d0; b1=1.085d-2; b2=-3.118d-5
       si0bc_ice=b0+(b1*T_ice)+(b2*T_ice*T_ice)-0.1d0
       del1bc_ice=cubicint_ice(Tc, T0bc, T0bc+5d0, 1d0, hbc)
       end if

       RETURN

       END SUBROUTINE prop_ice

!*************************************************************
!   Subroutine gauspdf (normalized width of the updraft distribution).
!************************************************************

       SUBROUTINE gausspdf(x, dp, sigmav_ice,miuv_ice,normv_ice)


       real*8, intent(in) :: x
       real*8 sigmav_ice,miuv_ice,normv_ice
       real*8, intent(out) :: dp

       sigmav_ice =max(sigmav_ice, 0.01)
       normv_ice =max(normv_ice, 0.01)


       dp=EXP(-0.5d0*(x-miuv_ice)*(x-miuv_ice)/sigmav_ice/sigmav_ice) /
     &sigmav_ice/sq2pi_par/(normv_ice + 0.001)


       RETURN


       END SUBROUTINE gausspdf


!*************************************************************
!   Function cubicint_ice (cubic interpolation between y1 and y2 within a and b).
!************************************************************

       real*8 function cubicint_ice(y, y1, y2, a, b)

       real*8, intent(in) :: y, y1, y2, a, b
       real*8 :: A_, B_, a0, a1, a2, a3, d, AUX

       if (y .le. y1) then
       d=a
       goto 5065
       end if

       if (y .ge. y2) then
       d=b
       goto 5065
       end if


       AUX=y2-y1
       A_=6d0*(a-b)/(AUX*AUX*AUX)
       B_=a+(A_*(y1*y1*y1)/6d0)-(A_*(y1*y1)*y2*0.5d0)

       a0=B_
       a1=A_*y1*y2
       a2=-A_*(y1+y2)*0.5d0
       a3=A_/3d0
       d=a0+(a1*y)+(a2*y*y)+(a3*y*y*y)


 5065    cubicint_ice=d


       end function cubicint_ice

!*************************************************************
!   Function dcubicint_ice (used in the PDA08 spectrum).
!************************************************************

       real*8 function dcubicint_ice(y, y1, y2, a, b)

       real*8, intent(in) :: y, y1, y2, a, b
       real*8 :: A_, a0, a1, a2, a3, d, AUX

       if (y .le. y1) then
       d=0
       goto 5065
       end if

       if (y .ge. y2) then
       d=0
       goto 5065
       end if


       AUX=y2-y1
       A_=6d0*(a-b)/(AUX*AUX*AUX)

       a1=A_*y1*y2
       a2=-A_*(y1+y2)*0.5d0
       a3=A_/3d0
       d=(a1)+(2d0*a2*y)+(3d0*a3*y*y)


 5065   dcubicint_ice=d


       end function dcubicint_ice

!*************************************************************
! Function PDG07 (simplified ice nucleation
!                     spectra according to Phillips et. al. 2007).
! si is supersaturation wrt ice and T is in K
!************************************************************

       real*8 function PDG07_ice(si, T)

       real*8, intent(in) :: si, T
       real*8 :: N

       if (T .le. 243d0)then
       N=1000d0*exp(-0.388d0)*(exp(3.88d0*si)-1d0)
       else
       N=60d0*exp(-0.639d0)*(exp(12.96d0*si)-1d0)
       end if

       PDG07_ice=N

       end function PDG07_ice




!*************************************************************
! Subroutine INSPEC_ice
!  Provides the Ice Nuclei concentration (m-3)
! and the chracteristic freezing threeshold, DSh (Barahona & Nenes 2009), at given
! si and T. The variable typeofspec_ice (integer) has the values
! 1 Meyers et. al. 1992
! 2  Phillips et. al. 2007
! 3  Barahona 2011
! 4  Phillips et. al. 2008 (simplifed)
! 5  Phillips et. al. 2013 (simplifed)
! si is supersaturation wrt ice and T is in K

!      Written by Donifan Barahona
!      donifan.o.barahona@nasa.gov

!************************************************************

       subroutine INSPEC_ice(six, N, Dsh,np_ice,norg_ice, sigorg_ice,
     &   dorg_ice, dbc_ice,sigbc_ice,
     &   kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
     & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
     & vpresw_ice,vpresi_ice)

       real*8, intent(in) :: six,np_ice,norg_ice, sigorg_ice,
     &  dorg_ice, dbc_ice,sigbc_ice,
     &  kdust_ice, kbc_ice, shdust_ice, shbc_ice,
     & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
     & si0bc_ice,nbc_ice,P_ice, T_ice
       real*8, intent(out) :: N, Dsh,Nhet_dep,Nhet_dhf,fdust_dep
       real*8 :: Nd, Nbc, aux, Si_, SW, del0, ddel0, fc, delw0, ddelw0,
     & SW0, Hdust, Hbc, Nbase, dNd, dNbc, dNbase, dH, dfc, Ndaux,
     & dNdaux, dNorg, Norg, Ndustaux, frac, aux2, Dx2, fdep, Ndep, Ndhf,
     & dNdep, dNdhf, si, dfrac, Ncdep_, Ncdhf_ , fglassy, Nglassy,
     & dNglassy, SIW, D_grid_bio, n_grid_bio,vpresw_ice,vpresi_ice

       real*8, dimension(3) :: sig_array, the_array, frac_array
       real*8, dimension(1:nbindust_ice) :: ndust_ice, sigdust_ice,
     &   ddust_ice

       real :: n_iw, DSh_s , nbc_s, dbc_s, Asolo
       real, dimension (nbindust_ice) :: ndust_s, ddust_s
       integer :: index, kindex, mode,nbindust_ice

       si=six
       Si_=si+1d0
       SW=Si_*vpresi_ice/vpresw_ice
       Nd = zero_par
       Nbc = zero_par
       Norg = zero_par
       Nglassy = zero_par

       dNd = zero_par
       dNbc = zero_par
       dNorg = zero_par
       dNglassy = zero_par

       if ((six .lt. 0.02) .or. (T_ice .gt. 270.0)) then
       N=0d0
       Dsh=si
       return
       end if

       if (SW .ge. 1.0) then
       SW=1.0
       Si_=vpresw_ice/vpresi_ice
       si=Si_-1.0
       end if
       SIW=vpresw_ice/vpresi_ice

       sig_array = 0.0
       the_array = 0.0
       frac_array = 0.0


       select case (typeofspec_ice)

       case(1)
       N=1000d0*exp(-0.639d0)*(exp(12.96d0*si)-1d0)
       Dsh=1d0/12.96d0

       case(2)
       N=PDG07_ice(si, T_ice)
       if (T_ice .le. 243d0)then
       Dsh=1d0/3.88d0
       else
       Dsh=1d0/12.96d0
       end if

       case(3)

       Ndustaux=0.0d0
       DO index=1, nbindust_ice

       Ndustaux=Ndustaux+ndust_ice(index)
       end do


       if (si .le.shdust_ice) then
       Nd=(si/shdust_ice)*Ndustaux*effdust_ice* exp(-kdust_ice*
     &(shdust_ice-si))
       dNd=Nd*((1d0/si)+kdust_ice)

       else
       Nd=Ndustaux*effdust_ice
       dNd=0d0
       end if


       if (si .le.shbc_ice) then
       Nbc=(si/shbc_ice)*nbc_ice*effbc_ice* exp(-kbc_ice*(shbc_ice-si))
       dNbc=Nbc*((1d0/si)+kbc_ice)
       else
       Nbc=nbc_ice*effbc_ice
       dNbc=0d0
       end if

       N=Nd+Nbc
       if (((dNd+dNbc) .gt. 0d0) .and. (N .gt. 0.0)) then
       Dsh=N/(dNd+dNbc)
       else
       Dsh=0.0
       end if

       case(4)



       SW0=0.97d0
       delw0=cubicint_ice(SW, SW0, 1d0, 0d0, 1d0)
       ddelw0=dcubicint_ice(SW, SW0, 1d0, 0d0, 1d0)

       Nbase=PDG07_ice(si, T_ice)


       if (T_ice .le. 243d0)then
       dNbase=3.88d0*Nbase
       else
       dNbase=12.96d0*Nbase
       end if


       del0=cubicint_ice(Si_, si0dust_ice, si0dust_ice+0.1d0, 0d0, 1d0)
       ddel0=dcubicint_ice(Si_, si0dust_ice, si0dust_ice+0.1d0, 0d0,
     & 1d0)

       fc=0.5d0*del1dust_ice*del0
       dfc=0.5d0*del1dust_ice*ddel0

       Hdust=fc+((1d0-fc)*delw0)
       dH=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))

       if (Hdust .gt. 1d0) then
       Hdust=1d0
       dH=0d0
       end if


       aux=(2d0/3d0)*Hdust*(Nbase/0.76d0)*pi_ice/5.0d-7/4d0
       aux2=(2d0/3d0)*pi_ice/0.76d0/5.0d-7/4d0

       Nd=0d0
       dNd=0d0

       DO index =1, nbindust_ice


       Dx2=ddust_ice(index)*ddust_ice(index)*ddust_ice(index)*0.52*
     &acorr_dust

       frac=0.5d0*(1d0-erfapp(-log(ddust_ice(index)/0.1e-6) /
     &sigdust_ice(index)/sq2_par))

       Ndaux=frac*ndust_ice(index)*(1d0-exp(-aux*Dx2))

       Nd=Nd+Ndaux
       Ndaux=(frac*ndust_ice(index)-Ndaux)
       dNdaux=Ndaux*((dH*Nbase)+(Hdust*dNbase))*aux2*Dx2

       dNd=dNd+dNdaux

       END DO




       del0=cubicint_ice(Si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)
       ddel0=dcubicint_ice(Si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)

       fc=0.5d0*del1bc_ice*del0
       Hbc=fc+((1d0-fc)*delw0)
       dfc=0.5d0*del1bc_ice*ddel0
       dH=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))


       if (Hbc .gt. 1d0) then
       Hbc=1d0
       dH=0d0
       end if

       frac=0.5d0*(1d0 -erfapp(-log(dbc_ice/0.1e-6) /sigbc_ice/
     &sq2_par))


       Dx2=dbc_ice*dbc_ice*dbc_ice*0.52*acorr_bc

       aux=((1d0/3d0)-0.06d0)*Hbc*(Nbase/0.76d0)*pi_ice/2.7d-7
       aux2=((1d0/3d0)-0.06d0)*pi_ice/0.76d0/2.7d-7

       Nbc=nbc_ice*frac*(1d0-exp(-aux*Dx2))
       dNbc=(nbc_ice*frac-Nbc)*((dH*Nbase)+(Hbc*dNbase))*aux2*Dx2




       frac=0.5d0*(1d0-erfapp(-log(dorg_ice/0.1e-6) /sigorg_ice/
     &sq2_par))

       Dx2=dorg_ice*dorg_ice

       aux=0.06d0*Hbc*(Nbase/0.76d0)*pi_ice/9.1d-7
       aux2=0.06d0*pi_ice/0.76d0/9.1d-7


       Norg=norg_ice*frac*(1d0-exp(-aux*Dx2))
       dNorg=(norg_ice*frac-Norg)*((dH*Nbase)+(Hbc*dNbase))*aux2*Dx2




       N=Nd+Nbc+Norg

       fdust_dep = Nd
       Nhet_dep = N


       if (.FALSE.) then

       if (T_ice .lt. 210.0) then
       Nglassy= min(0.01 +0.0045*(210.0 -T_ice), 0.1)
       fglassy= 7.7211*1e-3*Si_ - 9.2688*1e-3
       fglassy=min(fglassy, 3.3587e-3)
       fglassy=max(fglassy, 0.0)
       else
       Nglassy = 0.0
       fglassy = 0.0
       end if
       Nglassy = np_ice*Nglassy*fglassy
       dNglassy =Nglassy*7.7211*1e-3
       N=N+Nglassy

       end if



       if ((dNd+dNbc+dNorg+dNglassy) .gt. 0d0) then
       Dsh=N/(dNd+dNbc+dNorg+dNglassy)
       else
       Dsh=0.0
       end if


       case (5)



       D_grid_bio =dorg_ice



       n_grid_bio = 0.0

       if (is_gocart) then
       ndust_s = SNGL(frac_dust*ndust_ice)
       nbc_s = SNGL(frac_bc*nbc_ice)
       Asolo = SNGL(0.25d0*frac_org*norg_ice*pi_par*dorg_ice*dorg_ice)
       else


       DO index =1, nbindust_ice
       frac=0.5d0*(1d0-erfapp(log(0.1e-6/ddust_ice(index)) /
     &sigdust_ice(index)/sq2_par))
       ndust_s(index) = SNGL(frac*ndust_ice(index))
       end do

       frac=0.5d0*(1d0-erfapp(log(0.1e-6/dbc_ice) /sigbc_ice/sq2_par))

       nbc_s = SNGL(frac*nbc_ice)

       frac=0.5d0*(1d0-erfapp(log(0.1e-6/dorg_ice) /sigorg_ice/
     &sq2_par))*0.25d0

       Asolo = SNGL(frac*norg_ice*pi_par*dorg_ice*dorg_ice)
       end if

       ddust_s=SNGL(ddust_ice)
       dbc_s = SNGL(dbc_ice*1.0)

       call EMPIRICAL_PARAM_PHILLIPS(SNGL(Si_), SNGL(SIW), SNGL(SW), (/
     &ddust_s/), (/ndust_s/), 5, (/dbc_s/), (/nbc_s/), 1, (/
     &SNGL(D_grid_bio)/), (/SNGL(n_grid_bio)/), 1, Asolo, n_iw, DSh_s
     & ,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice)

       N=DBLE(n_iw)
       DSh=DBLE(DSh_s)

       case default
       N=zero_par
       Dsh=0.0
       end select

       if (Dsh .ge. si) then
       Dsh=si
       end if

       if (T_ice .gt. 255.0) then
       N=zero_par
       Dsh=zero_par
       end if

       end subroutine INSPEC_ice


!*************************************************************
!  Subroutine INimmersion
!  Provides the Immersion IN (for activated droplets) concentration at given T(K) according to Barahona et al. GMD, 2014.
!      Written by Donifan Barahona
!      donifan.o.barahona@nasa.gov
!========================================

       subroutine INimmersion(INconc, dINconc, wparcel,dbc_ice,sigbc_ice
     &  ,nbc_ice,fdust_imm,fdrop_dust,fdrop_bc,ndust_ice, sigdust_ice,
     &   ddust_ice,nbindust_ice,vpresw_ice,vpresi_ice,T_ice)

       real*8, intent (in) :: wparcel,dbc_ice,sigbc_ice,nbc_ice
     &  ,fdrop_dust,fdrop_bc
       real*8, intent (out) :: INconc, dINconc,fdust_imm
       real*8 :: Nd, Nd_unc, Nd_coa, Nbc, ahet, frac, dfrac, Naux,dNaux,
     & Tx, nssoot, nsdust, ninbkg, SX , dnsd, dnss, dNd, dNbc, coolr,
     & min_ns_dust, min_ns_soot, dninbkg,vpresw_ice,vpresi_ice,T_ice

       real*8, dimension(3) :: sig_array, the_array, frac_array
       real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
       integer :: index, kindex, nbindust_ice

       if (T_ice .lt. Thom) then
       INconc = zero_par
       dINconc = zero_par
       return
       end if

       if (T_ice .ge. To_ice) then
       INconc = zero_par
       dINconc = zero_par
       return
       end if

       Tx=T_ice -273.16
       coolr=5.0e-3*wparcel
       min_ns_dust= 3.75e6
       min_ns_soot= 3.75e9

       SX=(vpresw_ice/vpresi_ice)-1.0


       ninbkg = 0.0
       dninbkg=0.0





       nsdust= max(exp(-0.517*Tx + 8.934)-min_ns_dust, 0.0)
       dnsd = max(0.517*nsdust, 0.0)

       nssoot= max(1.0e4*exp(-0.0101*Tx*Tx - 0.8525*Tx + 0.7667)-
     &min_ns_soot, 0.0)
       dnss = max(-(-2.0*0.0101*Tx -0.8525)*nssoot, 0.0)

       Naux=zero_par
       dNaux=zero_par


       if (is_gocart) then

       DO index =1,nbindust_ice

       ahet= ahet_dust(index)
       Naux=(1.0-exp(-nsdust*ahet))*ndust_ice(index)+Naux
       dNaux = exp(-nsdust*ahet)*ndust_ice(index)*dnsd*coolr*ahet+dNaux
       END DO
       ahet = ahet_bc
       else

       DO index =1,nbindust_ice

       ahet= ddust_ice(index)*ddust_ice(index)*ddust_ice(index)*0.52*
     &acorr_dust* exp(4.5*sigdust_ice(index)*sigdust_ice(index))
       Naux=(1.0-exp(-nsdust*ahet))*ndust_ice(index)+Naux
       dNaux = exp(-nsdust*ahet)*ndust_ice(index)*dnsd*coolr*ahet+dNaux
       END DO

       ahet=dbc_ice*dbc_ice*dbc_ice*0.52*acorr_bc* exp(4.5*sigbc_ice*
     &sigbc_ice)
       end if


       Nd=Naux*fdrop_dust
       dNd=dNaux*fdrop_dust

       Nbc=nbc_ice*(1.0-exp(-nssoot*ahet))*fdrop_bc
       dNbc= nbc_ice*exp(-nssoot*ahet)*fdrop_bc*dnss*coolr*ahet




       INconc=Nbc+ Nd + ninbkg
       dINconc=dNbc+dNd + dninbkg
       fdust_imm = Nd


       end subroutine INimmersion

!
!=======================================================================================
!=======================================================================================
!=======================================================================================
!!====================================================================================
!            EMPIRICAL PARAMETERISATION (Phillips et al. 2013, JAS)
!            Code contributed by Vaughan Phillips, University of Leeds
! Implementation:   Donifan Barahona donifan.o.barahona@nasa.gov
!====================================================================================


      SUBROUTINE EMPIRICAL_PARAM_PHILLIPS(SI, SIW, SW, D_grid_dust,
     & n_grid_dust, ijstop_dust, D_grid_soot, n_grid_soot, ijstop_soot,
     & D_grid_bio, n_grid_bio, ijstop_bio, A_solo, n_iw, DSH,
     & Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice)
      implicit none
      real, intent(IN):: SI, SIW, SW, A_solo
      real*8, intent(IN):: P_ice, T_ice
      real, dimension(:), intent(IN):: D_grid_dust, n_grid_dust,
     & D_grid_soot, n_grid_soot, D_grid_bio, n_grid_bio
      integer, intent(IN):: ijstop_dust, ijstop_soot, ijstop_bio
!     integer ijstop_dust, ijstop_soot, ijstop_bio

      real :: nin_a_nuc_dust, nin_a_nuc_soot, nin_a_nuc_bio,
     & nin_a_nuc_solo, num_ic_dust_imm, num_ic_soot_imm, num_ic_bio_imm,
     & num_ic_solo_imm

      real, intent (inout) :: DSH, n_iw
      real*8, intent (out) :: Nhet_dep,Nhet_dhf,fdust_dep

      real :: dn_in_dust, dn_in_soot, dn_in_bio, dn_in_solo, dNall,
     & dNaux, naux, SS_w, dH_frac_dust, dH_frac_soot, dH_frac_solo, aux,
     & dfdep, temperature_K, P_SAT, ahet_aux



      REAL :: RHO_CFDC, BASE_DUST_OMEGA, BASE_SOOT_PHILIC_OMEGA,
     & BASE_BIO_OMEGA, ALPHA_DUST, ALPHA_SOOT, ALPHA_bio,
     & FRACTION_DEPNUCL_WARM_DUST, PIE, BASE_SOLO_OMEGA,
     & TEMP_MAX_DUST_DEGC, TEMP_MAX_SOOT_DEGC, TEMP_MAX_bio_DEGC,
     & GLASS_FRAC
      PARAMETER( BASE_DUST_OMEGA = 2.0e-6, BASE_SOOT_PHILIC_OMEGA = 1.e-
     &7, BASE_BIO_OMEGA = 0.89e-6, BASE_SOLO_OMEGA = 5.6e-5,
     & GLASS_FRAC = 0.5, ALPHA_DUST = 2./3., ALPHA_SOOT = 1./3. - 0.03,
     & ALPHA_bio = 0.03, RHO_CFDC = 50000./(287.*228.15),
     & FRACTION_DEPNUCL_WARM_DUST = 0.15, PIE = 3.1415926,
     & TEMP_MAX_DUST_DEGC = -10., TEMP_MAX_SOOT_DEGC = -15.,
     & TEMP_MAX_bio_DEGC = -2.)

      real :: FAC_CORRECT_RH = 2., rho_AIDA
      real:: H_frac_dust, n_in, n_in_dust, n_in_ultra, n_in_dust_ultra,
     & CIHENC_dust, ESW, ESI, SS_i, n_in_soot_ultra, H_frac_soot,
     & H_frac_bio, n_in_soot, n_in_bio, n_in_bio_ultra, CIHENC_soot,
     & CIHENC_bio, delta_Si, delta_T, delta_Sw, n_in_max, SS_iw, rho

      real :: H_frac_solO, RHI, n_in_solO, n_in_solO_star, CIHENC_solO,
     & Psi_solO

      real :: mu, S_i_0, RH_crit, S_i_w_warm, S_i_w_cold, S_i_w,
     & tc_HM_degC
      real :: S_w_0, dep_frac, n_in_hat, n_in_tilde
      real*4 :: dH1smooth
      real :: EPS = 0.622
      integer :: ij
!intrinsic :: exp, DEXP, SIZE, DBLE



!print *, SIZE(n_grid_dust(:))
      if(ijstop_dust .ne. SIZE(n_grid_dust)) stop 6366
      if(ijstop_soot .ne. SIZE(n_grid_soot)) stop 6366
      if(ijstop_bio .ne. SIZE(n_grid_bio)) stop 6366
!     ijstop_dust = SIZE(n_grid_dust)
!     ijstop_soot = SIZE(n_grid_soot)
!     ijstop_bio = SIZE(n_grid_bio)



!Naux=12.96 !default
      naux =0.0
!     Nall =dNaux
!     Sh=0.0
      n_iw=0.0
      nin_a_nuc_dust=0.0; nin_a_nuc_soot=0.0; nin_a_nuc_bio=
     &0.0; nin_a_nuc_solo=0.0
      num_ic_dust_imm=0.0; num_ic_soot_imm=0.0; num_ic_bio_imm=
     &0.0; num_ic_solo_imm=0.0
!     AnningC uncom following four lines
      n_in_dust=0.0; dn_in_dust=0.0;n_in_soot=0.0
      dn_in_soot=0.0; dn_in_bio=0.0; n_in_bio=0.0;
      dn_in_solO=0.0;n_in_solO=0.0;n_in_dust_ultra=0.0
      n_in_soot_ultra=0.0;n_in_bio_ultra=0.0;


      H_frac_dust = 0.0
      H_frac_soot = 0.0
      H_frac_solo = 0.0
!     H1smooth=0.0
      aux=0.0

      temperature_K=SNGL(T_ice)
      P_SAT =SNGL(P_ice)

!A_solo = 1e-7 !m2 kg-1


!====================================================================================
!            COMPUTATION BLOCK
!
!====================================================================================
!
       rho_AIDA = 90000./(287.*205.)

       rho = P_SAT/(287.*temperature_K)

       Psi_solO = A_solO/BASE_SOLO_OMEGA
       SS_i = min(max(SI-1.0, 0.0), 1.0)
       SS_w = min(max(SW-1.0, -1.0), 1.0)
       SS_iw = min(max(SIW - 1.0, 0.0), 1.0)


       if(SS_i > 0.0) then
       if(temperature_K < 273.15 .and. temperature_K > 273.15 -
     & 90. ) then


       if(SS_w > 0.) then
       SS_i = SS_iw
       SS_w = 0.0
       end if


!                        S_i_zero = 1.15 !this is taken care of


       delta_Si = H_1_smooth(SS_i + 1, 1.1, 1.2, 0.0, 1.,dH1smooth);
       delta_T = H_1_smooth(-(temperature_K-273.15), 35., 40.,
     & FRACTION_DEPNUCL_WARM_DUST, 1.,dH1smooth);
       delta_Sw = H_1_smooth(SS_w + 1.0, 0.97, 1., 0., 1.,dH1smooth);

       tc_HM_degC = temperature_K - 273.15


       S_i_0 = 1. + 10.**(8.2584e-6*tc_HM_degC*tc_HM_degC*tc_HM_degC +
     & 5.3938E-4*tc_HM_degC*tc_HM_degC + 3.1656E-3*tc_HM_degC - 1.0261)


       S_w_0 = 0.97

       aux =H_1_smooth(-(temperature_K-273.15), 35., 40.,
     & FRACTION_DEPNUCL_WARM_DUST, 1.,dH1smooth)/FAC_CORRECT_RH

       dep_frac = H_1_smooth(SS_i + 1, S_i_0, S_i_0 + 0.1, 0.,1.,
     &  dH1smooth)* aux
       dfdep=dH1smooth*aux

       aux= H_1_smooth(SS_w + 1.0, S_w_0, 1., 0.,1.,dH1smooth)

       H_frac_dust = dep_frac + (1. - dep_frac)*aux

       dH_frac_dust = dfdep + (SIW*(1. - dep_frac)*dH1smooth)- aux*
     &dfdep

       if(H_frac_dust > 1.) H_frac_dust = 1.

       if ((H_frac_dust .gt. 1.0e-6) .and. (H_frac_dust .lt. 1.)) then
       dH_frac_dust = dH_frac_dust/H_frac_dust
       else
       dH_frac_dust =0.0
       end if


       S_i_0 = 1.2

       aux =H_1_smooth(-(temperature_K-273.15), 65., 75., 0.,1.,
     &   dH1smooth)
       dep_frac = H_1_smooth(SS_i + 1, S_i_0, S_i_0+0.1, 0.,1.,
     &    dH1smooth)*aux
       H_frac_solO = dep_frac
       dH_frac_solo=0.0
       if ((H_frac_solo .gt. 1.0e-6) .and. (H_frac_solo .lt. 1.)) then
       dH_frac_solo = dH1smooth/H_frac_solo
       end if
       if(H_frac_solO > 1.) H_frac_solO = 1.


       S_w_0 = 0.97

       S_i_0 = 1.3
!
       aux = H_1_smooth(-(temperature_K-273.15), 40., 50., 0.,1.,
     &   dH1smooth) / FAC_CORRECT_RH
       dep_frac = H_1_smooth(SS_i + 1, S_i_0, S_i_0+0.1, 0.,1.,
     &   dH1smooth)* aux

       dfdep= dH1smooth*aux

       aux = H_1_smooth(SS_w + 1.0, S_w_0, 1., 0.,1.,dH1smooth)
       H_frac_soot = dep_frac + (1. - dep_frac)*aux
       if(H_frac_soot > 1.) H_frac_soot = 1.

       dH_frac_soot = dfdep + (SIW*(1. - dep_frac)*dH1smooth)- aux*
     &dfdep
       if ((H_frac_soot .gt. 1.0e-6) .and. (H_frac_soot .lt. 1.)) then
       dH_frac_soot = dH_frac_soot/H_frac_soot
       else
       dH_frac_soot =0.0
       end if


       H_frac_bio = H_frac_soot

       if(temperature_K < 273.15 .and. temperature_K >= 273.15 -
     & 35.) then
       n_in = 1.E3* (exp(12.96*SS_i - 0.639)/RHO_CFDC) *0.0587*
     &FAC_CORRECT_RH
       if( temperature_K > 273.15 -5. .and. temperature_K < 273.15 -
     & 2. ) then
       n_in = n_in*H_1_smooth(-(temperature_K-273.15), 2., 5., 0., 1.,
     &   dH1smooth)
       endif
       if(temperature_K >= 273.15 - 2. ) n_in = 0.


       if(temperature_K < 273.15 -25. ) then
       n_in_tilde = 1000.*(exp(0.1296*(SS_i*100.-10.))**0.3)*
     &FAC_CORRECT_RH/RHO_CFDC
       n_in_hat = n_in

       if(temperature_K >= 273.15 - 30.) n_in_max = 1.E3* (exp(12.96*
     &SS_iw - 0.639)/RHO_CFDC) *0.0587*FAC_CORRECT_RH
       if(temperature_K < 273.15 - 30.) n_in_max = 1000.*(exp(0.1296*
     &(SS_iw*100.-10.))**0.3)*FAC_CORRECT_RH/RHO_CFDC

       if(n_in_hat > n_in_max) n_in_hat = n_in_max
       if(n_in_tilde > n_in_max) n_in_tilde = n_in_max



       n_in = n_in_hat * ((n_in_tilde/n_in_hat)**(H_1_smooth(-
     &(temperature_K-273.15), 25., 35., 0., 1.,dH1smooth)))


       if(n_in > n_in_max) n_in = n_in_max

       endif
       n_in_dust = 0.
       dn_in_dust = 0.


       if(temperature_K < 273.15 - 30.) then
       dnaux = 3.88
       else
       dnaux = 12.96
       end if


       naux=0.0
       do ij = 1, ijstop_dust
       if (is_gocart) then
       mu = n_in*ALPHA_DUST*H_frac_dust*ahet_dust(ij)/BASE_DUST_OMEGA
       else
       ahet_aux = PIE*D_grid_dust(ij)*D_grid_dust(ij)*D_grid_dust(ij)*
     &4.73*acorr_dust/6.0
       mu = n_in*ALPHA_DUST*H_frac_dust*ahet_aux/BASE_DUST_OMEGA
       end if
       naux = (1. - exp(-mu))*n_grid_dust(ij)
       n_in_dust = n_in_dust + naux
       dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +
     & dH_frac_dust), 0.0) + dn_in_dust

       enddo

       if( temperature_K > 273.15 +TEMP_MAX_DUST_DEGC -
     & 20. .and. temperature_K < 273.15 + TEMP_MAX_DUST_DEGC) then
       n_in_dust = n_in_dust*H_1_smooth(-(temperature_K-273.15),-
     &TEMP_MAX_DUST_DEGC,-TEMP_MAX_DUST_DEGC+20., 0., 1.,dH1smooth)
       endif
       if(temperature_K >= 273.15 + TEMP_MAX_DUST_DEGC) n_in_dust = 0.


       n_in_soot = 0.
       dn_in_soot = 0.
       do ij = 1, ijstop_soot

       if (is_gocart) then
       mu = n_in*ALPHA_SOOT*H_frac_soot*ahet_bc/BASE_SOOT_PHILIC_OMEGA
       else
       ahet_aux = PIE*D_grid_soot(ij)*D_grid_soot(ij)*D_grid_soot(ij)*
     &16.4*acorr_bc/6.0
       mu = n_in*ALPHA_SOOT*H_frac_soot*ahet_aux/BASE_SOOT_PHILIC_OMEGA
       end if
       naux = (1. - exp(-mu))*n_grid_soot(ij)
       n_in_soot = n_in_soot + naux
       dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux+dH_frac_soot),
     & 0.0) + dn_in_soot

       enddo

       if( temperature_K > 273.15 + TEMP_MAX_SOOT_DEGC -
     & 10. .and. temperature_K < 273.15 + TEMP_MAX_SOOT_DEGC) then
       n_in_soot = n_in_soot*H_1_smooth(-(temperature_K-273.15),-
     &TEMP_MAX_SOOT_DEGC,-TEMP_MAX_SOOT_DEGC+10., 0., 1.,dH1smooth)

       endif
       if(temperature_K >= 273.15 + TEMP_MAX_SOOT_DEGC) n_in_soot = 0.

       n_in_bio = 0.
       dn_in_bio = 0.


       do ij = 1, ijstop_bio
       mu = n_in*ALPHA_bio*H_frac_bio*PIE*(D_grid_bio(ij)**2.) /
     &BASE_BIO_OMEGA

       mu = n_in*ALPHA_bio*H_frac_bio
       naux = (1. - exp(-mu))*n_grid_bio(ij)




       n_in_bio = n_in_bio + naux
       dn_in_bio = max(mu*(n_grid_bio(ij)-naux)*dnaux, 0.0) + dn_in_bio


       enddo


       if( temperature_K > 273.15 + TEMP_MAX_bio_DEGC -
     & 3. .and. temperature_K < 273.15 + TEMP_MAX_bio_DEGC) then
       n_in_bio = n_in_bio*H_1_smooth(-(temperature_K-273.15),-
     &TEMP_MAX_bio_DEGC,-TEMP_MAX_bio_DEGC+3., 0., 1.,dH1smooth)

       endif
       if(temperature_K >= 273.15 + TEMP_MAX_bio_DEGC ) n_in_bio = 0.



       else
       n_in = 0.; n_in_ultra = 0.; n_in_dust = 0.; n_in_soot =
     & 0.; n_in_bio = 0.;
       endif

       if(temperature_K < 273.15 - 35.) then
       n_in_ultra = 1000.*(exp(0.1296*(SS_i*100.-10.))**0.3)*
     &FAC_CORRECT_RH/RHO_CFDC
       dnaux = 3.88
       naux=0.0


       RHI = (SS_i+1.)*100.
       if(RHI < 0.) RHI = 0.
       n_in_solO_star = 1000.e6*(7.7211e-5 * RHI - 9.2688e-3)/rho_AIDA


       n_in_dust_ultra = 0.;
       dn_in_dust = 0.0
       do ij = 1, ijstop_dust

       if (is_gocart) then
       mu = n_in_ultra*ALPHA_DUST*H_frac_dust*ahet_dust(ij)/
     &BASE_DUST_OMEGA
       else
       ahet_aux = PIE*D_grid_dust(ij)*D_grid_dust(ij)*D_grid_dust(ij)*
     &4.73*acorr_dust/6.0
       mu = n_in_ultra*ALPHA_DUST*H_frac_dust*ahet_aux/BASE_DUST_OMEGA
       end if



       naux = (1. - exp(-mu))*n_grid_dust(ij)
       n_in_dust_ultra = n_in_dust_ultra + naux
       dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +dH_frac_dust),
     & 0.0) + dn_in_dust



       enddo



       n_in_soot_ultra = 0.0
       dn_in_soot = 0.0
       do ij = 1, ijstop_soot

       if (is_gocart) then
       mu = n_in_ultra*ALPHA_SOOT*H_frac_soot*ahet_bc/
     &BASE_SOOT_PHILIC_OMEGA
       else
       ahet_aux = PIE*D_grid_soot(ij)*D_grid_soot(ij)*D_grid_soot(ij)*
     &16.4*acorr_bc/6.0
       mu = n_in_ultra*ALPHA_SOOT*H_frac_soot*ahet_aux/
     &BASE_SOOT_PHILIC_OMEGA
       end if



       naux = (1. - exp(-mu))*n_grid_soot(ij)
       n_in_soot_ultra = n_in_soot_ultra + naux
       dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux +dH_frac_soot),
     & 0.0) + dn_in_soot

       enddo


       n_in_bio_ultra = 0.
       dn_in_bio = 0.0










       n_in_solO = Psi_solO*glass_frac*H_frac_solO*n_in_solO_star
       dn_in_solO =max(Psi_solO*glass_frac* (H_frac_solO*77211.0*100.0/
     &rho_AIDA + n_in_solO_star*dH_frac_solo), 0.0)


       else
       n_in_ultra = 0.; n_in_dust_ultra = 0.; n_in_soot_ultra =
     & 0.; n_in_solO = 0.; n_in_bio_ultra = 0.;
       endif



       n_in_dust = n_in_dust + n_in_dust_ultra;
       n_in_soot = n_in_soot + n_in_soot_ultra;
       n_in_bio = n_in_bio + n_in_bio_ultra;




! PROBLEM:  how to ensure that the frozen fraction does not exceed 1 ?

       if (.false.) then
       if(n_in_dust + n_in_bio + n_in_soot + n_in_solO > 0.) then

       CIHENC_dust = n_in_dust - nin_a_nuc_dust
       if(CIHENC_dust < 0.) CIHENC_dust = 0.

       CIHENC_soot = n_in_soot - nin_a_nuc_soot
       if(CIHENC_soot < 0.) CIHENC_soot = 0.

       CIHENC_bio = n_in_bio - nin_a_nuc_bio
       if(CIHENC_bio < 0.) CIHENC_bio = 0.

       CIHENC_solO = n_in_solO - nin_a_nuc_solO
       if(CIHENC_solO < 0.) CIHENC_solO = 0.


       n_iw = n_iw + CIHENC_dust
       nin_a_nuc_dust = nin_a_nuc_dust + CIHENC_dust
       num_ic_dust_imm = num_ic_dust_imm + CIHENC_dust

       n_iw = n_iw + CIHENC_soot
       nin_a_nuc_soot = nin_a_nuc_soot + CIHENC_soot
       num_ic_soot_imm = num_ic_soot_imm + CIHENC_soot

       n_iw = n_iw + CIHENC_bio
       nin_a_nuc_bio = nin_a_nuc_bio + CIHENC_bio
       num_ic_bio_imm = num_ic_bio_imm + CIHENC_bio

       n_iw = n_iw + CIHENC_solO
       nin_a_nuc_solO = nin_a_nuc_solO + CIHENC_solO
       num_ic_solO_imm = num_ic_solO_imm + CIHENC_solO


       endif
       end if
       endif
       endif

       n_iw = n_in_dust + n_in_bio + n_in_soot + n_in_solO
       dNall = dn_in_dust + dn_in_bio + dn_in_soot + dn_in_solO

       if (n_iw .gt. zero_par) then
       fdust_dep = DBLE(n_in_dust/n_iw)
       end if
       Nhet_dep = DBLE(n_in_dust)
       Nhet_dhf = DBLE(n_in_solO)




       if (( dNall > 0.) .and. (n_iw .gt. 0.0)) then
       Dsh=max(min(n_iw/dNall, SS_i), 0.005)
       else
       Dsh=0.005
       end if


      END SUBROUTINE EMPIRICAL_PARAM_PHILLIPS

      real function H_1(X, X_1, X_2, Hlo)
      real, intent(in) :: Hlo, X, X_1, X_2

      if(X >= X_2) H_1 = 1
      if(X <= X_1) H_1 = Hlo
      if(X > X_1 .and. X < X_2) H_1 = (X - X_1)/(X_2 - X_1)

      if( X_2 <= X_1) stop 91919

      return
      end function


      real function H_1_smooth(X, X_1, X_2, Hlo, Hhi,dH1smooth)
      real, intent(in) :: Hlo, Hhi, X, X_1, X_2
      real*4,intent(out) :: dH1smooth
      real :: a_0, a_1, a_2, a_3, A, B

      if(X >= X_2) H_1_smooth = Hhi
      if(X <= X_1) H_1_smooth = Hlo

      if(X >= X_2) dH1smooth = 0.0
      if(X <= X_1) dH1smooth = 0.0

      if(X > X_1 .and. X < X_2) then
       A = 6.*(Hlo - Hhi)/(X_2**3. - X_1**3. + 3.*(X_2*X_1*X_1 - X_1*
     &X_2*X_2) )
       a_3 = (A/3.)
       a_2 = -(A/2.)*(X_1 + X_2)
       a_1 = A*X_2*X_1
       B = Hlo + A*(X_1**3.)/6. - A*X_1*X_1*X_2/2.
       a_0 = B
       H_1_smooth = a_0 + a_1*X + a_2*X*X + a_3*X*X*X
       dH1smooth = a_1 + 2.0*a_2*X + 3.0*a_3*X*X
      endif

!H1smooth =min(dH1smooth , 1000000.0)
!H1smooth =max(dH1smooth , 0.0)

      if( X_2 <= X_1) stop 91919

      return
      end function





!    END ICE PARAMETERIZATION DONIF
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC



       END MODULE aer_cloud
